6576
Langmuir 2007, 23, 6576-6587
Brownian Dynamics Simulations of Associating Diblock Copolymers M. J. Cass and D. M. Heyes* DiVision of Chemistry, School of Biomedical and Molecular Sciences, UniVersity of Surrey, Guildford GU2 7XH, United Kingdom
R. J. English North East Wales Institute, UniVersity of Wales, Plas Coch Campus, Mold Road, Wrexham LL11 2AW, Wales, United Kingdom ReceiVed NoVember 2, 2006. In Final Form: March 20, 2007 A novel coarse-grained computational model for associating polymers is proposed that is based on a Gaussian “blob” representation of the polymer chains. The model allows a large number of model polymers to be simulated at moderate computational cost over a wide packing fraction range using the Brownian dynamics, BD, technique. The attraction of the hydrophobic part of the polymer to those on other molecules can lead to strong aggregation of the polymer molecules in real systems, and this is included in the model by an attractive potential felt by the Gaussian blobs to a common “nodal” point that represents the center of the micelle. Attention here is confined to model AB diblock copolymers in which the hydrophilic block, A, has a much higher mass than the hydrophobic moiety, B, which leads to relatively small aggregation numbers, Nagg, of ∼8. The aggregation number at low packing fractions is found to increase with packing fraction, as observed in experiments, with a functional form that closely follows a simple theory derived here that is based on entropy-derived mean-field terms for the free-energy change associated with the incorporation of the polymer molecule into the micelle. The computational model exhibits an extremely low critical micelle concentration (cmc), and micelles with Nagg ≈ 5 are observed at the lowest packing fractions, φ, simulated (∼10-4), which is consistent with experiment. The long-time self-diffusion coefficient of the polymers (and hence micelles) decreases logarithmically with packing fraction, and the viscosity increased with concentration according to the Huggins equation. The spherical blob coarse graining results in the simulable time scales being longer than the Rouse time of the chain, and hence for the nonassociating polymers the intrinsic viscosity is an input parameter in the model. The introduction of association leads to the partial inclusion of the intrinsic viscosity in the simulation and has an effect on the computed Huggins coefficient, kH, which is found to be ∼6 in those cases.
I. Introduction The field of complex liquids has attracted much attention not only because of the rich phase behavior but also because of the significant commercial and biological importance of these systems. Within this field, associating polymers, APs, constitute an important class. These polymers contain “sticker” groups that are less soluble in the solvent and cause association. Generally the polymer consists of hydrophilic and hydrophobic blocks. The hydrophobic blocks act as sticker groups that can order the polymers into supramolecular structures, even at low concentrations. (We will refer exclusively here to water-based systems, but the same principles could be extended to other solvents.) APs are widely employed in industry, primarily as rheology modifiers,1 but they also have uses elsewhere, ranging from drug delivery and other bioapplications2-5 to templating nanoparticles.6 The rheology of an associating polymer solution depends on chemical composition and molecular architecture, with the number of sticker groups in the polymer having perhaps the greatest * Corresponding author. E-mail:
[email protected]. (1) Polymers As Rheology Modifiers; Schultz, D. N., Glass, J. E., Eds.; Advances in Chemistry Series 462; American Chemical Society: Washington, DC; Oxford University Press: Oxford, U.K., 1991. (2) Alexandridis, P.; Lindman, B. Amphiphilic Block Copolymers: SelfAssembly and Applications; Elsevier: Amsterdam, 2000. (3) Wang, W.; McConaghy, A. M.; Tetley, L.; Uchegbu, I. F. Langmuir 2000, 17, 631. (4) Jiang, G.-B.; Quan, D.; Liao, K.; Wang, H. Mol. Pharm. 2006, 2, 152. (5) Malmsten, M. Surfactants and Polymers in Drug DeliVery; Taylor & Francis: Oxford, U.K., 2002. (6) Hamley, I. W. Introduction to Soft Matter John Wiley & Sons: Chichester, U.K., 2000.
affect. Single sticker APs can form micellar structures similar to those produced by surfactants. Spherical micelles form when the size of the hydrophobe is much less than that of the hydrophile. With two stickers (telechelics), the self-assembly behavior has more possibilities because individual polymers can have both ends in the same micelle (loop arrangement) or two different micelles (bridge state), allowing the potential for the formation of networks of micelles and percolating gel-like structures. With many stickers in the molecule, more complex structures including intramolecular micelles can also form.7 For linear polymers, it is typical for the sticker to be located at one end of the polymer (a diblock copolymer, AxBy) or in the center (a triblock copolymer, AxByAx, with the hydrophobic group in the middle), for example, where the hydrophilic group, A, is poly(oxyethylene) and the hydrophobic group, B, is poly(oxypropylene). The aggregation number and critical micelle concentration (cmc) of (EO)x(PO)y(EO)x depend on the polymer molecular mass, the ratio x/y, and the temperature.8 The temperature dependence of the cmc is derived from that of the hydrophobicity of the PO block.8-10 The cmc can be extremely low; for example, PStPEp in n-decane at 30° has a cmc of 2.4 × 10-8g cm-3,11 and for (EO)x(PO)y(EO)x, values as low as 10-5 (7) Borisov, O. V.; Halperin, A. Curr. Opin. Colloid Interface Sci. 1998, 3, 415. (8) Wanka, G.; Hoffmann, H.; Ulbricht, W. Macromolecules 1994, 27, 4145. (9) Alexandridis, P.; Holzwarth, J. F.; Hatton, T. A. Macromolecules 1995, 27, 2414. (10) Jebari, M. M.; Ghaouar, N.; Aschi, A.; Gharbi, A. Polym. Int. 2006, 55, 176. (11) Price, C. Pure Appl. Chem. 1983, 55, 1563.
10.1021/la063210j CCC: $37.00 © 2007 American Chemical Society Published on Web 05/12/2007
Brownian Dynamics Simulations of Diblock Copolymers
% concentration by mass were found.9 The self-assembly of block copolymers into micelles is similar to that of surfactants, although surfactant micelle behavior is dominated by packing considerations whereas the entropy of the polymer chain is more important for the coblock polymers and scaling theory would be expected to apply. There is a significant body of literature on scaling theories of the phase behavior and rheological properties of APs (e.g., refs 12 and 13), which is based on mean-field representations of the key processes and time scales in the system. These treatments are successful in reproducing the observed experimental trends14-16 because they depend on chain statistics on the Kuhn length scale and longer. Molecular computer simulation is a complementary approach that is capable of representing the self-assembly of these systems at the molecular level. It is able to provide insights into the nature of self-assembly and various static properties and dynamical and rheological properties as well as relationships between them. The approach has fewer assumptions that the purely analytic theories because it explicitly includes manybody effects and is therefore quite distinct in its starting point from scaling theories. A fully atomistic simulation treatment is currently not feasible because of the large number of atoms in these systems and the long time scales associated with the network dynamics and rheology (compared to characteristic atomistic motion time scales). The interactions between a large number of polymer chains must be simulated, and with polymers aggregating into micelles, many such micelles need to be incorporated in the computer model to minimize system size effects. To reduce the computational costs, a coarse-grained or mesoscale model for the diblock molecule is used. There have been a number of coarse-grained models developed to simulate block copolymer systems on length and time scales much larger than accessible using atomistic models. (See ref 17 for a review of these methods.) These involve representing many atoms by a single sphere and the block copolymer by a small number of these spheres or “beads” held together by springlike interactions along the chain. These models have been able to reproduce not only spontaneous micellization but also the various microdomain phases exhibited by these systems over a range of concentrations and temperatures.18-20 We have also developed mesoscale models for polymers in which the polymer molecule is represented by a small number of Kuhn beads,21 in a similar vein to Groot and Agterof’s Monte Carlo studies of associative polymer networks22,23 in which the AP molecule was represented by a small number of hard beads. By comparison, the model presented here uses a much softer effective polymer-polymer interaction potential without any specific number of beads. It employs a free-energy interaction potential between the moieties and allows for multiply associating stickers. As stated earlier, the simplest single-sticker AP architectures are of the diblock copolymer type with A and B equal-length (12) Rubinstein, M.; Semenov, A. N. Macromolecules 2001, 1058, 34. (13) Annable, T.; Buscall, R.; Ettelaie, R. Colloids Surf., A 1996, 112, 97. (14) Annable, T.; Buscall, R.; Ettelaie, R.; Whittlestone, D. J. Rheol. 1993, 37, 695. (15) Molino, F.; Appell, J.; Filali, M.; Michel, E.; Porte, G.; Mora, S.; Sunyer, E. J. Phys: Condens. Matter 2000, 12, A491. (16) Pellens, L.; Ahn, K. H.; Lee, S. J.; Mewis, J. J. Non-Newtonian Fluid Mech. 2004, 121, 87. (17) Ortiz, V.; Nielsen, S. O.; Klein, M. L.; Discher, D. E. J. Polym. Sci., Part B 2006 44, 1907. (18) Groot, R. D.; Madden, T. J. J. Chem. Phys. 1998, 108, 8713. (19) Groot, R. D.; Madden, T. J.; Tildesley, D. J. J. Chem. Phys. 1999, 110, 9739. (20) Liu, W.; Qian, H.; Lu, Z.; Li Z.; Sun, C. Phys. ReV. E 2006, 74, 021802. (21) Xiao, C.; Heyes, D. M. J. Chem. Phys. 1999, 111, 10694. (22) Groot, R. D.; Agterof, W. G. M. J. Chem. Phys. 1994, 100, 1649. (23) Groot, R. D.; Agterof, W. G. M. J. Chem. Phys. 1994, 100, 1657.
Langmuir, Vol. 23, No. 12, 2007 6577
blocks or where one block is much larger than the other. These two extremes are more conveniently simulated than architectures in between. In either case, coarse graining is most appropriate for high-molecular-weight molecules. The properties of micelles of a polymer with roughly equal-length A and B blocks would be expected to be dominated by the packing of the hydrophobic, B, blocks in the micelle core. If the hydrophilic group, A, is much larger than B, then the polymers form “hairy” micelles with comparatively small aggregation numbers,24,25 which is advantageous for simulation because a large number of micelles can then be considered. Additionally, such asymmetric block polymers are useful prototypes for multisticker AP architectures such as telechelics and comb polymers in which the hydrophobic groups are also generally much smaller than the hydrophilic backbone.26-28 The advantage of coarse graining on the scale of the radius of gyration of the polymer (which is assumed to be equivalent to that of the hydrophilic moiety) is that this makes it feasible to follow the dynamical evolution of the polymer network and its rheology for sufficiently long time scales to capture the important dynamical processes taking place in these systems and obtain satisfactory statistical averaging.
II. Theory A. Computational Model. The first issue to be considered is the force field between polymers. Because the model is based on a mesoscale representation of the diblock copolymers, these interaction laws should be thought of as free energy rather than enthalpic pair interactions. Excluded Volume. In this study, attention is focused on model AB diblock copolymers in which the (insoluble) hydrophobe (B) is much smaller than the hydrophilic moiety (A). The AP is represented by an excluded volume, EV, center of interaction coupled to an attractive group at a point in space. The core-core interactions have a single Gaussian analytic form, as developed by Louis and Hansen,29-32 which is the most coarse-grained representation of a polymer molecule. Support for the single blob model for a highly asymmetric diblock is provided by smallangle neutron-scattering experiments on swollen block copolymer micelles, where the interaction between the micelles can be represented by the hard sphere potential.33 The rheology of block copolymer micelles in the range of 10-100 nm has also been modeled as systems of soft spheres.34 In fact, the representation of polymers of various architectures in solution by soft interaction laws is a well-established preocedure in theoretical models for these systems.35 The characteristic length scale of the model is much larger than the Kuhn length and is taken to be the radius of gyration, (24) Stellbrink, J.; Rother, G.; Laurati, M.; Lund, R.; Willner, L.; Richter, D. J. Phys: Condens. Matter 2004, 16, s3821. (25) Forster, S.; Zisenis, M.; Wenz, E.; Antonietti, M. J. Chem. Phys 1996, 104, 9956. (26) Hamley, I. W. Block Copolymers in Solution; John Wiley & Sons: Chichester, U.K., 2005. (27) Tirtaatmadja, V.; Tam, K. C.; Jenkins, R. D. Macromolecules 1997, 30, 3271. (28) Tam, K. C.; Seng, W. P.; Jenkins, R. D.; Bassett, D. R. J. Polym. Sci., Part B 2000, 38, 2019. (29) Louis, A. A.; Bolhuis, P. G.; Hansen, J. P. Phys. ReV. E 2000, 62, 7961. (30) Lang, A.; Likos, C. N.; Watzlawek, M.; Lowen, H. J. Phys: Condens. Matter 2000, 12, 5087. (31) Louis, A. A.; Bolhuis, P. G.; Hansen, J. P.; Meijer, E. J. Phys. ReV. Lett. 2000, 85, 2522. (32) Hansen, J. P.; Addison, C. I.; Louis, A. A. J. Phys: Condens. Matter 2005, 17, S3185. (33) Castelletto, V.; Hamley, I. W.; English, R. J.; Mingvanish, W. Langmuir 2003, 19, 3229. (34) Buitenhuis, J.; Fo¨rster, S. J. Chem. Phys 1997, 107, 262. (35) Likos, C. N.; Hoffmann, N.; Lo¨wen; Louis, A. A. J. Phys.: Condens. Matter 2002, 14, 7681.
6578 Langmuir, Vol. 23, No. 12, 2007
Cass et al.
rgyr, of the polymer. The number of Kuhn segments in the polymer is not explicitly specified in the model but is assumed to be large enough that the polymer is in the asymptotic scaling regime. (Typically, they have about 100 Kuhn lengths.) The polymer is treated as a single blob. Such a model has proved to be quite successful at concentrations up to ca. c*, the overlap concentration,36,37 but it clearly cannot account for entanglement effects and reptation that dominate the system’s physical properties at higher concentrations. The EV potential between the AP cores is
UEV(r) )
{
EEV exp(-r2/2σEV2) r < 3σEV r > 3σEV 0
(1)
Figure 1. Possible schemes for introducing hydrophobic association into a computer model. The regions in bold indicate the location of the hydrophobic interaction. (a) Coarse graining of the associating polymers with distinct interaction centers for the hydrophobic and hydrophilic moieties. (b) Coarse graining of the entire polymer as a single central interaction point from which the potential energy has a minimum at intermediate range. (c) Model employed in this study, where the associative polymer ends are introduced only when attached to a node.
where EEV is the magnitude of the excluded volume interaction at maximum overlap (i.e., r ) 0), which is approximately 2kBT, kB is Boltzmann’s constant, and T is the temperature.29 The characteristic length σEV is equal to the radius of gyration, rgyr, of the polymer.29 Scaling theories for micelles with a small core and large corona predict that the aggregation number and core radius are independent of coronal chain length,26 which gives support for this level of coarse graining in which there is no formal inclusion of the number of Kuhn length beads. Equation 1 can be viewed as a potential of mean force, an idea explored by others.38 If the polymer is subdivided into more segments, the multibead and spring model is produced, which is more typical of previous simulations of polymers. The intra- and intermolecular beadbead interactions could still be represented by the potential form of eq 1 but with increasing barrier height (i.e., UEV(0)) because the chain structure is more finely resolved. Such a model was employed by Xiao and Heyes in Brownian dynamics simulations of unassociating polymers in various quality solvents21 and for associating polymers39 in which EEV ≈ 50kBT was used for a 12-bead model for the polymer. Such a highly coarse-grained representation of the polymer molecule as a single or a very small number of Gaussian blobs serves as the basis for the present model of associating polymers. For a symmetric diblock or comparable moiety triblock polymer, it would be reasonable to represent each block as a distinct Gaussian blob, with the mutual interactions between the blobs being given by eq 1. Hansen et al.32,40 also showed that diblock copolymers with equal block sizes can be represented by two repulsive Gaussian beads joined by a spring potential. It is important to note that representing each polymer block as a separate Gaussian blob works well when the blocks are similar sizes but can present computational difficulties when there is a significant size difference between blocks, as will be discussed in the next section. In addition, the fact that we use a spherically averaged Gaussian potential for the interaction between the chains implies that the time step in the simulation is larger than the Rouse relaxation time of the isolated chain; therefore, the model cannot give information on processes taking place on sub-Rouse time scales (e.g., for the intrinsic viscosity, see below). Hydrophobic Interactions. The incorporation of the hydrophobe in the interaction model posed the greatest challenge in the model construction, and various types of hydrophobe-hydrophobe and hydrophobe-hydrophile interactions were explored. Some pos-
sible ways of incorporating the hydrophobic interactions are illustrated in Figure 1. In the first example shown in Figure 1a, a large, purely repulsive blob representing the hydrophilic part, A, of the molecule is attached by a “spring” potential to a small blob representing the hydrophobic block, B. These B groups attract each other with a potential that can give rise to association of the molecules. The molecule has two units executing BD equations of motion, with a much smaller friction coefficient associated with the latter because its hydrodynamic radius is smaller. The time step has to be very small to track the hydrophobes (compared to the optimum value for the hydrophile), and too many time steps are required to explore phase space adequately. An alternative strategy that was investigated was to treat the hydrophobe as a modification of the EV potential representing the polymer. A single Gaussian core was surrounded by an attractive “well” at intermediate separation (in the spirit of the Lennard-Jones potential41). This approach, represented schematically in Figure 1b, failed to form micellar structures. The attraction between polymers was insufficiently localized because the attractive potential was spherically symmetric about the polymer center of mass rather than being directional in nature. In addition, the range of possible analytic forms for attraction was restricted to ensure compliance with Ruelle’s most necessary but not sufficient condition for thermodynamic stability (that the integral of the potential over all space must be positive).29,42,43 If Ruelle’s condition is not obeyed, then there is insufficient penalty for overlap of the model polymer molecules and the total energy of the system would continually decrease as more polymer molecules become localized into a small region of space. The low critical micelle concentration observed for the AP molecules suggests that the hydrophobes spend the vast majority of their time associated with other hydrophobes.44 The adopted model utilizes this observation. The hydrophobes are only explicitly included in the model when they are associated with other hydrophobes. The model is illustrated in Figure 1c. The AP molecules are attracted toward a central point or “node” that represents the average position of the center of the micelle. The node is created when at least two AP molecules come within a specified range. Additional molecules can become attached to the node to increase the size of the micelle. The node is removed from the system when there are fewer than two molecules
(36) Bolhuis, P. G.; Louis, A. A.; Hansen, J. P.; Meijer, E. J. J. Chem. Phys. 2001, 144, 4296. (37) Bolhuis, P. G.; Louis, A. A. Macromolecules 2002, 35, 1860. (38) Akkermans, R. L. C.; Briels, W. J. J. Chem. Phys. 2001, 114, 1020. (39) Xiao, C.; Heyes, D. M. J. Chem. Phys. 2002, 117, 2377. (40) Addison, C. I.; Hansen, J. P.; Krakoviack, V.; Louis, A. A. Mol. Phys. 2005, 103, 3045.
(41) Heyes, D. M. The Liquid State; John Wiley & Sons: Chichester, U.K., 1998; p 64. (42) Ruelle, D. Statistical Mechanics: Rigorous Results; Imperial College Press: London, 1969. (43) Fisher, M. E.; Ruelle, D. J. Math. Phys. 1966, 7, 260. (44) Tripathi, A.; Tam, K. C.; McKinley, G. H. Macromolecules 2006, 39, 1981.
Brownian Dynamics Simulations of Diblock Copolymers
Langmuir, Vol. 23, No. 12, 2007 6579
Figure 2. Key behavior in the evolution of nodes. (a) Two unassociated polymers are close to each other (r < 2ra), and their hydrophobes associate to form a node. (b) An unassociated polymer becomes attached to an existing node when its separation from the node is less that ra, and associated beads will dissociate from nodes when their separation exceeds ra. (c) Two nodes are in close proximity (r < ra/2) and merge to form a larger micelle. rf ) 2ra is the “formation” distance between two EV cores when they form a node.
associated with it. The hydrophobe position is unspecified unless it is attached to a node, when it is assumed to be located at the node. Because the relaxation time scale of the hydrophobes and, by implication, of the nodes is much faster than that of the EV units, it is appropriate to evolve them in a different manner. The position of a node represents the average location of the associated hydrophobes; therefore, its position at each time step was set to the average position of the EV cores. This location should nevertheless generally be close to the minimum potential energy position for the node because all of the interaction potentials in the model are slowly varying with distance. The BD algorithm was applied only to the EV units, thereby permitting a relatively large time step. As implied above, when two unassociated polymers came close enough together to form a node, the position of the node was assumed to be midway between them. For more than two APs, the node was located at the center of mass of the bound polymer molecules. The polymers were attracted to the node via a potential, Ua(r), where r is the separation of the center of mass of the AP from the node (i.e., micelle core). Ua(r) represents the combined effects of hydrophobe-hydrophobe attraction and the extension of the hydrophobic end of the polymer from its center of mass. An attractive Gaussian analytic form was chosen for Ua(r) (the precise potential form is not known so a Gaussian is a reasonable choice)
Ua(r) )
{
-Ea exp(-Rr2) r < ra r > ra 0
(2)
where Ea is the free energy of association and ra is the range of association/dissociation from the node. The value of R was set to 9/2ra2 so that the potential falls to a negligible value when r ) ra (ra is then 3 standard deviations away from the center of the Gaussian). ra could be taken as the approximate radius of the core of the micelle. At small r, the polymer experiences a nearquadratic or harmonic potential field as is expected from entropic considerations. Harmonic spring models of polymers, such as the FENE potential45 at small extension, are widely used in polymer chain modeling at the coarse-grained level. At larger separation, the potential falls to zero, which represents on this scale of description the gradual dissociation or detachment of (45) Zhou, Q.; Akhavan, R. J. Non-Newtonian Fluid Mech. 2004, 116, 269.
the polymer from the node. A frequently used free-energy model of micelle formation suggests that Ea should vary with micelle size, m (m is the number of hydrophobes associated with a node). The average surface area of a hydrophobe in contact with the water phase is dependent on the size of the micelle core and hence m.46 For surfactants, the Gibbs energy of association per molecule, ∆G0m/m ≈ Ea, is usually taken to have the form46,47
∆G0m ) -B1 + B2m-1/3 + B3m1/3 m
(3)
where B1, B2, and B3 are positive constants. B1 is the free-energy gain associated with the removal of the hydrophobe from solution into the micelle, and B2 is the residual surface energy per hydrophobe in contact with the water in the micelle. B3 represents the repulsion between the coronal components of the amphiphilies. The last term does not need to be included in the node-polymer potential because this term represents the excluded volume repulsion of polymer molecules in the micelle so it is already explicitly included in the computer model through the potential in eq 1. Therefore, the association potential in the current model is equivalent to
∆G0m ) -B1 + B2m-1/3 m
(4)
where ∆G0m/m can be equated approximately to Ea. Trial simulations in which Ea was varied according to eq 4 gave a micelle size distribution that was hardly different from those for which Ea was assumed to be independent of m. Hence to maintain the simplicity of the model, Ea was treated as a variable parameter (see the discussion below) but was not treated as a function of the micelle size. As mentioned earlier, the association and dissociation of polymers in the model is characterized by the depth of the node potential, Ea, and its range, ra. There are three distinct processes that need to be considered, as illustrated in Figure 2. These are the formation of a node from two unassociated polymers, association with (and dissociation from) a pre-existing node, and (46) Warr, G. G.; White, L. R. J. Chem. Soc., Faraday Trans. 2 1985, 81, 549. (47) Leibler, L.; Orland, H.; Wheeler, J. C. J. Chem. Phys. 1983, 79, 3550.
6580 Langmuir, Vol. 23, No. 12, 2007
the merging of two nodes. An unassociated polymer is assumed to associate with a node when its separation from that node is less than ra. Polymers associated with a node are taken to dissociate from that node when their separation from it exceeds ra, with the potential between the node and bead having fallen to essentially zero at this point. If the number of polymers associated with a node falls below two, then the node ceases to exist. The formation of nodes occurs when the COM of two unassociated polymers come within 2ra of each other, forming a node midway between them. The resulting node is thus less than ra from either polymer. If two nodes come close enough to each other, then it is physically realistic for them to merge to form a single node because the hydrophobes in an equivalent real system would be drawn together. It is important that the range at which nodes merge is significantly less than ra to avoid unphysical dissociation of polymers from merging nodes. (The separation of the center of mass of a polymer could be greater than ra from the position of the new node.) By trial and error, the range for the merging of nodes was taken to be ra/2. It is possible for an unassociated polymer COM to be within the range of several nodes after a move, in which case a choice between the possible nodes for association must be made. This selection was weighted by the relative reduction in the polymer’s potential energy for the various possible associations. The probability of a particular association was proportional to the potential energy reduction of that association. Of two possible associations if one had an attractive potential energy gain twice the other, that association was twice as likely. The values of ra and Ea for eq 2 are the only free parameters in the model, and they were selected within a physically reasonable range by trial and error. The effects of the range of choices are discussed in the Results section. It is worth highlighting what the model does not include. As a result of the coarse graining, much of the chemical detail relevant to the atomistic-scale behavior is not retained, although one’s expectation is that those chemical features are retained that have an influence on the mesoscale behavior. There is no explicit water, and the model polymers are treated as moving in a continuous background. The representation of the micellar core does not have any oligomeric detail. No allowance is made for the interfacial free energy between the core and the water. Details of the chemical linkage between the blocks are not included. The core is represented by a node to which the polymer molecules are attracted via a central force interaction that decays with distance. Therefore, the sticker has a variable sticking energy depending on the distance of the center of mass of the polymer from the node. The excluded volume interaction between the chains is assumed to be spherically symmetric and concentrationindependent. The simple representation of the core means that we cannot directly fit our model within the usual packing parameter description used for oligomeric surfactant systems. No allowance is made for chain deformation in the micelle as is predicted in, for example, the analytic Daoud-Cotton model.48 Louis et al. investigated compaction effects on the effective potential between the centers of mass of polymers in solution.29 They calculated the effective excluded volume potential between polymers as a potential of mean force by inverting the radial distribution function between the centers of mass from Monte Carlo simulations of self-avoiding walk (SAW) polymer chains on a lattice. The height of the potential barrier increased by only about 10% from infinite dilution up to about 4 times the overlap concentration (approximately the range covered in the present study). The range of the interaction did not change significantly (48) Daoud, M.; Cotton, J. P. J. Phys. (Paris) 1982, 43, 531.
Cass et al.
in this concentration regime either. In the micelle, there is overlap and stretching of the chains. The precise analytic form for the excluded volume interactions is still not known in the micellar geometry. We have assumed a state point independent form for this interaction for this pilot study. The computations were performed with a free-draining model for the hydrodynamics. Although many-body hydrodynamics (MBH) effects could be included,49 they add considerably to the computer time, and this would limit the number of polymers that could be simulated. It is possible that MBH are not very important for the associative polymers because the physics is dominated by the thermodynamics of self-assembly. In fact, as will be shown, the present simple model does manifest the essential features and trends seen in experimental studies of asymmetric coblock polymers. The main purpose of this work is to test the node construction model for the self-assembly of associative polymers, which is a new methodology. Refinements to the model to include (a) distortion of the core, (b) chain stretching of the corona and more complicated excluded volume interactions, and (c) manybody hydrodynamics, although of course highly desirable as the next stage, are not necessary (even appropriate) for this preliminary study. B. Technical Details and Computed Quantities. The system was evolved using a standard free-draining Brownian dynamics, BD, algorithm,50,51 which gives both static and dynamic properties. Each time step the EV cores were moved by BD, and then the node positions were recalculated and associations were made or broken as appropriate. The simulations were carried out in reduced units of the radius of gyration, rgyr ) 1, β-1 ) kBT ) 1, and solvent viscosity, ηs ) 1, so that in these units EEV ) 2 and σEV ) 1. The reduced time step used was ∆t ) 0.5, which was established from ∆t ) β〈∆r2〉3πηsrh, where 〈∆r2〉 is the meansquare displacement in the time step and rh is the hydrodynamic radius that was taken to be 0.8rgyr, a value appropriate to that of a good solvent.52 This gives a mean-square displacement for each polymer of ∼0.07 per time step. Simulations were carried out for a range of packing fractions in the range of φ ) 10-4 to 1 for typically 30 000 time steps in each case. The polymer Newtonian zero shear viscosity, ηp, is the key parameter for the AP systems, and only the Newtonian viscosity is considered in this study. ηp was calculated from the off-diagonal components of the pressure tensor, pRβ,
pRβ )
1 V
rijR fijβ ∑i ∑ j>i
(5)
where fijβ is the β Cartesian component of the force between particles i and j using a Green-Kubo equation
η - ηs ) ηp ) βV
∫0∞ 〈pRβ(s) pRβ(s + t)〉s dt
(6)
where ηs is the solvent viscosity, V is the volume of the simulation cell, and 〈‚‚‚〉s indicates an average over time origins, s. Only the excluded volume and node forces are included in eq 5 because all other interactions are lost in the coarse-graining procedure and the solvent is included only indirectly via the Brownian motion term. As a result of this coarse graining, intrapolymer interactions are present only when polymers are associated through the potential, Ua, given in eq 2, which includes the attraction (49) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352. (50) Bran´ka, A. C.; Heyes, D. M. Phys. ReV. E 1998, 58, 2611. (51) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: San Diego, CA, 2002. (52) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press: Oxford, U.K., 2003.
Brownian Dynamics Simulations of Diblock Copolymers
Langmuir, Vol. 23, No. 12, 2007 6581
between the chain end and its center of mass. The intrinsic viscosity, [η], gives the contribution to the viscosity in the limit of infinite dilution from the polymer’s internal degrees of freedom. For the unassociated polymers, the model does not generate any contribution to the viscosity from the intrinsic viscosity, and even when polymers are associated, the intrinsic viscosity is only partially included. This is a consequence of the coarsegraining procedure in which a spherical blob is adopted, so the average shape of the polymer is taken as the starting point and hence the relevant time scales are beyond the Rouse time (which characterizes the time scale of shape fluctuations of the polymer). Consequently, the calculated viscosity reflects mainly the second and higher terms in the usual virial or Huggins expansion of the polymer solution viscosity
η - ηs ) [η] + kH[η]2c + ‚‚‚ ηsc
(7)
where c is the concentration of polymer and kH is the Huggins coefficient,52 which is dimensionless. The first term on the right side of eq 7, [η], represents the contribution to the viscosity from the internal degrees of freedom of the polymer, and the second term involving kH reflects the intermolecular interactions. In the absence of an associative interaction, eq 6 gives only the intermolecular contribution to the viscosity, ηip, rather than ηp (the polymer’s explicit contribution to the viscosity). The viscosity of the real system can be written, however, as the sum of these two terms
η - ηs ) [η]c + ηip ηs
(8)
η - ηs ) ηp ηs
(9)
rather than
(Note that, as defined, ηip is the intermolecular contribution to the viscosity divided by the solvent viscosity and is hence dimensionless.) To compare the simulation results with experiment, a transferable concentration scale needs to be established, which can be achieved by benchmarking against the experimental viscosity of real nonassociating polymers. Comparing eq 8 with the Einstein equation for the shear viscosity of monodispersed spheres, η/ηs ) 1 + 2.5φ + ‚‚‚,53 where φ is the packing fraction, it can be seen that 2.5φ ) [η]c. Thus from eqs 7 and 8,
kH )
ηip (2.5φ)2
(10)
The above relationship can be used to relate a convenient definition of the packing fraction for the simulation data to the experimental value. A reasonable definition of the effective packing fraction of the polymers, φ′, is
φ′ )
Nπ(2 ln 2)3/2 6V
(11)
where N is the number of polymers simulated and (2 ln 2)1/2 is the value of r (in units of rgyr) when the excluded volume potential of eq 1 is equal to 1kBT. This is the separation of maximum force between the blobs. By use of eq 10, the ratio of the simulated (53) van de Ven, T. G. M. Colloidal Hydrodynamics; Academic Press: London, 1989; p 222.
packing fraction to the experimental value can be found. The computed viscosity in the simulation depends, to leading order, on the square of φ′ because there is no computed intrinsic viscosity contribution (for nonassociating polymers). Hence,
η - ηs ) ηip ) k′H(φ′)2 ηs
(12)
where k′H is the Huggins coefficient found from the simulation data. A comparison of eqs 12 and 10 using 5φ/2 ) c[η] gives k′H ) ηip/(φ′)2. Hence, k′Hφ′2 ) 6.25kHφ2; therefore,
x
2 φ ) φ′ 5
k′H kH
(13)
where kH is the experimental Huggins coefficient, which for nonassociating polymers in a good solvent is in the relatively narrow range of 0.3 < kH < 0.4.52 This typically gives φ ≈ 2φ′. The formula in eq 13 uses the relationship between the low density behavior of the computed viscosity and that obtained experimentally.
III. Results and Discussion Real-time visualizations of the polymer system from the computer model confirmed that micellization occurred in broad agreement with expectations. Micellar structures were produced, and the network evolved significantly on a typical simulation time scale. Polymers could be seen to associate and dissociate from nodes. Figure 3 presents a few representative snapshots from a simulation. These snapshots are taken from simulations with a small number of polymers because with a large number it is difficult to distinguish individual micelles on a 2D projection of a 3D simulation box. The straight lines in the Figure represent the links between a polymer EV and a node through the potential defined in eq 2. Figure 3a is at a low packing fraction of φ′ ) 5 × 10-3. Figure 3b is at the higher packing fraction of φ′ ) 5 × 10-2, and Figure 3c is a different perspective of the configuration of part b. Although the micelles are on average spherical, instantaneous configurations show significant deviations from spherical symmetry. In any simulation configuration, there is a distribution of micelle sizes measured by m, the number of polymers in a given micelle. The time-averaged micelle size distribution, P(m), is shown in Figure 4a. In Figure 4b, the average aggregation number, Nagg, ∞ ∞ m2P(m))/(∑m)1 mP(m)), is shown. In both where Nagg ) (∑m)1 a and b, these are shown as a function of system size for the same packing fraction (φ′ ) 2 × 10-2). The probability distribution and Nagg are quite sensitive to N for small N. The right-hand frame shows the time-dependent buildup of the aggregation number. (Each simulation started from a random arrangement of model polymer molecules.) This Figure illustrates the importance of simulating a sufficiently large number of polymers. With only 100 polymers, Nagg at long times was ∼14, whereas for N ) 2000 the long-time limiting value is ∼12. The distribution of micelle sizes can also be seen to be skewed to larger values when an insufficient number of polymers are simulated. The micelle size distributions were statistically the same for simulations of 500, 1000, and 2000 polymers. A thousand model polymers were used in the simulations for the data presented in the Figures below, which provided results with minimal system size artifacts without unduly increasing the simulation running time.
6582 Langmuir, Vol. 23, No. 12, 2007
Cass et al.
Figure 3. Examples of the micelles formed in simulation. The small gray patches represent the individual polymer molecules, and the lines denote the links between polymers and the micelle core or node. Key: (a) low packing fraction (φ′ ) 5 × 10-3); (b and c) two views of the same higher packing fraction simulation (φ′ ) 5 × 10-2). For clarity, these pictures were produced using only 50 model polymers in the simulation, whereas 1000 polymers were used in the production simulations to obtain the system averages used for most of the remaining Figures.
As stated earlier, the simulation packing fraction must be scaled to match the experimentally equivalent value, φ ) 0.4[η]c, through the Huggins coefficient,
k′H )
ηip ηs(φ′)2
(14)
Figure 5 shows that from the simulation ηip values k′H ) 7.5 ( 1.0 in the limit of zero density. By use of eq 13, the relationship between the experimental and simulation system packing fractions is φ ) (1.85 ( 0.27)φ′ on taking kH ) 0.35 (0.3 < kH < 0.4 for polymers in a good solvent52). Figure 6 shows the dependence of the micelle size probability distribution, P(m), on the association/dissociation distance, ra, used in eq 2 for Ea ) 14kBT. As the value of ra increases from 4.5 to 6.5, the distribution broadens and shifts to larger m values. A simple mean-field statistical mechanical, MF, theory for micelle formation and P(m) can be derived that uses an approximate analytic expression for the energy of a polymer molecule in the micelle, U(m). This includes contributions from the excluded volume repulsion between polymers forming the micelle and attractive polymer-
Figure 4. Comparison of (a) the micelle size, m, distribution, P(m), and (b) the average micelle size, Nagg () 〈m〉) for simulations with different numbers of model polymer molecules in the simulation cell, N. Key: N ) 100 (-), N ) 500 (‚‚‚), N ) 1000 (---), and N ) 2000 (-‚-). These simulations were all carried out for φ′ ) 2 × 10-2 with Ea ) 14kBT and ra ) 5.5.
node interactions. The average excluded volume energy of polymer i interacting with the other polymers associated with the node, UEV(i), is
〈 ∑ ( )〉 -rij2
m-1
UEV(i) )
EEV exp
j*i
2
(15)
in reduced distance units using eq 1. rij is the separation between the model polymer molecules i and j, and 〈‚‚‚〉 indicates as usual a statistical average. Assume that the EV cores are randomly
Brownian Dynamics Simulations of Diblock Copolymers
Langmuir, Vol. 23, No. 12, 2007 6583
Figure 5. Concentration-dependent estimate of the Huggins coefficient using k′H ) ηip/(φ′)2 for the case of nonassociating coarsegrained polymers (i.e., with Ea ) 0) using eq 14. The Huggins coefficient is this ratio in the limit of zero concentration (---). Error bars are shown for each state point. The error bars in k′H are derived from the difference between the viscosities of two independent simulations.
Figure 7. Same as for Figure 6 except that the dependence of the probability distribution on Ea is shown for a constant value of the association distance, ra ) 5.5. The distributions from left to right are for Ea ) 10kBT, 14kBT, and 18kBT. The corresponding curves from the MF approach using the same parameters are also shown.
The total energy, U(i), of the EV core i including the attractive term with respect to the micelle center from eq 2 is
U(i) )
EEV(m - 1) (1 - exp(-2rs2)) - Ea exp(-Rrs2) (18) rs2
In the model, the energy of a polymer in a micelle of size m is U(m) ≡ U(i), which is a function of EEV, Ea, m, ra, and rs. For a given set of the first four parameters, rs can be found by numerically minimizing the total energy of the micelle, ∂U(i)/∂rs ) 0. The probability function, P(m), is given by
P(m) )
Figure 6. Micelle size probability distribution P(m), where m is the number of hydrophobes in a micelle as function of the range of association (ra) for Ea ) 14kBT. From left to right, the distributions are shown for ra ) 4.5, 5.5, and 6.5. An intermediate value of the packing fraction, φ′ ) 2 × 10-2, was used in each case, with 400 polymers simulated. Also shown are the corresponding distributions predicted by the mean-field, MF, approach using the same parameters as discussed in the text associated with eq 19.
distributed over a spherical surface of radius rs about the node. The distance between two randomly chosen points on this surface is 2r(1 - cos(θ)), where θ is the angle between the lines to the two points from the micelle center of mass. Performing an integration over the surface using standard spherical coordinates gives
UEV(i) ) EEV
∫0∞∫0π∫02π F(r) 2 exp(-r2(1 - cos(θ)))r2 sin(θ) dφ dθ dr
) 4πEEV
∫0∞ F(r)(1 - exp(-2r2)) dr
(16)
where the number density of EVs is F(r) at r. The EV cores are assumed to be randomly distributed on a spherical shell at a fixed distance rs from the node, hence F(r) ) δ(r - rs)(m 1)/4πr2 and
UEV(i) )
EEV(m - 1) (1 - exp(-2rs2)) 2 rs
(17)
exp(-βU(m))
∑n
(19) exp(-βU(n))
To simplify the treatment, it is assumed that every monomer has the same energy in the system (i.e., the system is monodisperse in micelles of size m, so there is no degeneracy term included in eq 19). This procedure can be used to obtain P(m) for given values of ra and Ea. As can be seen in Figure 6, this gives a very similar trend for the dependence of P(m) on ra to that observed from the simulation. As would be expected, given that rs is fixed to a single value in the MF treatment, the distributions seen in simulation are broader, and this effect increases with ra. The larger ra, the greater the variation in the node-polymer separation and the less valid the assumption that all polymers are at the same separation from the node. Figure 7 shows the dependence of the micelle size probability distribution function, P(m), on the node well depth, Ea, of the potential in eq 2. Ea values in the range (10-18)kBT were used in the simulations. The distribution P(m) shifts to higher m values with increasing Ea, whereas the peak height hardly changes. As for Figure 6, the MF model and simulations manifest the same qualitative trends. Figure 8 shows the radial density of the EV cores about the nodes as a function of Ea for simulations at φ′ ) 2 × 10-2. The peak position in these profiles decreases and the magnitude increases a little with increasing strength of attraction. However, the Figure suggests that there is a relatively small dependence on the size of the micelle with Ea/kBT in the range of 10-18. The changes in the average separation of polymers from their node with Ea are consistent with MF theory, although a larger Ea leads to a greater distribution of m in the simulation, which is incompatible with the MF model.
6584 Langmuir, Vol. 23, No. 12, 2007
Figure 8. Node-polymer probability density distribution from simulation as a function of Ea. The plots show increasing Ea from 10kBT to 18kBT in steps of 2kBT. The maximum probability increases with Ea and moves to smaller values of the node-core distance, rnb. The simulation parameters are φ′ ) 2 × 10-2, N ) 1000, and ra ) 5.5.
Cass et al.
F127 ((EO)106(PO)70(EO)106), aggregation numbers of about 10 were recorded close to the critical micelle temperature (m ) 10 for PE6800 at 30 °C and m ) 7 for F127 at 20 °C). These simulations are based on the assumption that the EV core has spherical symmetry. One would expect the real chains to adopt a more prolate ellipsoidal shape with increasing m to minimize their overlap, an effect that should increase the average value of m. The critical micelle concentration of the model appears to be very low, and significant aggregation is evident in Figure 9 even at the low density of φ ) 1.85 × 10-4, where φ is the dimensionless packing fraction (φ ) 2/5c[η]). Leibler et al.47 developed a mean-field theory that includes the concentration dependence of the micellar aggregation number. Although in their case this dependence was weak, it is quite possible that for the smaller micelles generated in these simulations the φ dependence is more significant. To investigate this further, a simplified version of their theory is developed below. The key to establishing the concentration dependence is the inclusion of the translational entropy, SM, of a “gas” of micelles in the expression for the free energy, FM, of the micellar phase.
FM ) NmF - TSM
(20)
where FM is the free energy of the micellar phase, Nm is the number of micelles of size m, F is the free energy of a micelle, and T is the temperature. To simplify the problem, the expression of Warr and White for F from ref 46 is used, and it is assumed that the contribution to the free energy from the interface between the micellar core and water is negligible (i.e., B2 ) 0), giving
F ) -B1m + B3m4/3 - (B3 - B1)
Figure 9. Average micelle size, Nagg, as a function of packing fraction for the simulation parameters N ) 1000, Ea ) 16, and ra ) 6.0. φ ) 1.85φ′. For each time step in the simulation, the average micelle size in the system was recorded, and the error bars on the graph were derived from the minimum and maximum values of these during the simulation.
Figure 9 shows the average micelle size, Nagg ) 〈m〉, as a function of packing fraction for the parameter set N ) 1000, Ea ) 16, and ra ) 6. The micelle size increases monotonically with packing fraction over several orders of magnitude, and there is a steeper rise near φ ≈ 0.1. The Figure shows that there are essentially two regions. At low concentration, there is a roughly logarithmic dependence of Nagg on concentration up to φ ≈ 0.01. The average micelle sizes are very similar to those found in Monte Carlo simulations of diblock copolymers employing a multibead representation of the individual polymer molecules by Timoshenko.54 The aggregation numbers from the simulation are at the low end of what is observed experimentally for diblock systems, which is possibly because the micelle core in the current model is assumed to have negligible volume. The aggregation number found experimentally varies considerably between different chemical systems, but for those cases with a large x/y ratio, the values are relatively small and similar to those found here. For example, Fo¨rster et al.25 found aggregation numbers of 6 and 14 for their samples PS - 2,1 and PS - 3,1 respectively. For triblock pluronics such as PE6800 ((EO)74(PO)30(EO)74) and (54) Timoshenko, E. G.; Kuznetsov, Y. A. Europhys. Lett. 2001, 53, 322.
(21)
where the constants Bi are the same as those defined for eq 3. The final term in eq 21 ensures that F ) 0 for m ) 1. The micelle size is obtained by minimizing the total free energy with respect to m,
dFM dNmF dSm )0) -T dm dm dm
(22)
The translational entropy is given by
SM 1 ) -Nm ln(φζ) + - 1 ln(1 - φζ) kB φζ
(
(
)
)
(23)
where ζ is the fraction of polymers in micelles. Taking ζ ≈ 1 and using Nm ) N/m (which assumes that all micelles are the same size, i.e., m ) Nagg), where N is the number of polymers in the system, gives
dSM N 1 ) kB 2 ln(φ) + - 1 ln(1 - φ) dm φ m
(
(
)
)
(24)
Thus, using eq 24 with eq 21 gives
(
)
(B3 - B1) dNmF d )N -B1 + B3m1/3 dm dm m )N
(
)
B3 -2/3 (B3 - B1) m + 3 m2
(25)
Brownian Dynamics Simulations of Diblock Copolymers
Langmuir, Vol. 23, No. 12, 2007 6585
Figure 10. Plot of Nagg4/3 against φ on a logarithmic scale. The square symbols with the error bars are the simulation Nagg, and the dashed line is from eq 28. The least-squares fit parameters are B1 ) 13.0, and B3 ) 1.18, and φ ) 1.85φ′. The diamonds (temperature of 60 °C) and other square symbols (temperature of 20 °C) are data taken from Borbe´ly, S. Langmuir 2000, 16, 5540 on the aggregation number obtained by small-angle neutron scattering on Brij-35 (PEObased nonionic surfactants). The concentration scale for these experimental data is on the top axis. The error bars were assigned as for Figure 9.
Noting that kBT ≡ 1 in reduced units,
N
(
)
B3 -2/3 (B3 - B1) + ) m 3 m2
(
)
N 1 ln(φ) + - 1 ln(1 - φ) (26) φ m2
(
)
Hence,
m4/3 )
3B1 3 1 ln(φ) + - 1 ln(1 - φ) - 3 + (27) B3 φ B3
(
(
)
)
For low concentrations, the approximation ((1/φ) - 1) ln(1 - φ) ) -1 can be made, to give
Nagg4/3 ) m4/3 ≈
( )
3 (ln(φ) - 1 - B3 + B1) B3
that
(28)
If micelles are to form, then the magnitude of the free-energy loss on removal of a hydrophobe from water must be greater than the gain associated with the repulsion between polymers in the micelle. The greater the water-removal free energy, the larger the micelles and the lower the concentration at which they will form. If B1 - B3 - 1 is small (say 10 rather than 100), then Nagg will be sensitive to φ, whereas if this constant is much larger, Nagg will be relatively insensitive to φ. For these model micelles, we have relatively small aggregation numbers, suggesting that B1 - B3 - 1 e10. The theory predicts that m4/3 is a linear function of ln(φ) for low φ, which because this simple treatment assumes monodisperse micelles (i.e., Nagg ) m) predicts that Nagg4/3 is a linear function of ln(φ) for low φ, which is tested in Figure 10. Figure 10 gives Nagg4/3 against φ for the parameters N ) 1000, Ea ) 16.0, and ra ) 6.0. The initial slow increase in aggregation number at low concentrations, as evident in Figure 9, is consistent with eq 28. The density variation of micelle size is sensitive to the loss of entropy when micelles are formed. In systems with a predisposition to large aggregation numbers (i.e., B1 . B3), eq
Figure 11. Potential per monomer, u, against packing fraction. The simulation parameters are N ) 1000, Ea ) 16, and ra ) 6.0. φ ) 1.85φ′.
28 indicates a weak dependence of Nagg on concentration.26,47,55 However, with the small aggregation numbers (i.e., B1 ≈ 3B3 + 3) found here the concentration dependence is more significant. A fit to the simulated data can be seen in Figure 10, giving values of B1 ) 13.0 and B3 ) 1.18 that are consistent with this explanation. The fitted values are also compatible with the simulation model parameters. It would be expected that B1 ≈ Ua(〈rnb〉), where rnb is the node-bead separation, and B3 ≈ UEV(〈rbb〉) where rbb is the separation of nearest-neighbor polymer beads within the micelle corona. Assuming that rnb ≈ ra/3 ) 2.0 and rbb ≈ σEV ) 1.0, then B1 ) 16 exp(-0.5) ) 9.7 and B3 ) 2.0 exp(-0.5) ) 1.2, which are quite close to the values extracted from the linear plot in Figure 9. Also shown in Figure 10 are the aggregation number data obtained by small-angle neutron scattering on Brij-35 (PEO-based nonionic surfactants). It can be seen that these data also follow the Nagg4/3 against φ predicted trend quite well. Quantitative agreement between experiment and simulation could presumably be achieved by fine tuning the model parameters. The steep rise in aggregation number at high packing fraction, as evident in Figure 9, occurs where a significant overlap of the polymer chains would occur in real systems. The average potential energy per polymer, 〈u〉, is positive in this packing fraction range, as shown in Figure 11, m
〈u〉 )
{Ua(i) + ∑UEV(rij)} ∑ i)1 j*1 N
(29)
where Ua(i) is the node potential experienced by polymer i, UEV(i, j) is the excluded volume potential between polymers i and j, and N is the number of polymers in the system. When the potential per polymer becomes positive, there is a significant degree of overlap of the Gaussian cores, and the present node construction model is most severely tested. Further into the semidilute regime, the correlation length, ξ, decreases, and each polymer would be better represented by a necklace of blobs, the size of which increasingly tends to the Kuhn segment length. The self-diffusion coefficient of the monomers, D, was computed from their mean-square displacement, 〈(δr)2〉, using D ) (∂〈(δr)2〉/∂t)/6 in the limit of t f ∞.51 The same equation was used to compute the self-diffusion coefficient of the micelle, Dm, where the mean-square displacement of the node was (55) Brinke, G. T.; Hadziioannou, G. Macromolecules 1987, 20, 486.
6586 Langmuir, Vol. 23, No. 12, 2007
Figure 12. Plot of D/D0 against packing fraction φ, where D0 ) 1.0/(6πx3-(2/Nagg)) in reduced units; we have assumed that the hydrodynamic radius is equal to the radius of gyration for an ideal star polymer with Nagg arms and that a monodisperse micelle distribution can be used. The simulation parameters are N ) 1000, Ea ) 16, and ra ) 6.0. φ ) 1.85φ′. The line is meant to guide the eye. The error bars are calculated from the distribution of the slopes of the mean-square displacement of the model polymer center of mass against time (for D) and from those in D0. The standard error in the value of D0 () 1/6πηsrh), where rh is the hydrodynamic radius of the micelle, was taken from that of the average aggregation number of the micelle, as for Figure 9.
considered instead. The statistics of the node are poorer than for the polymers because there are more of the latter and the finite node lifetime limits the number of samples that can be taken for the micelle statistics. It was confirmed that the value of D for the micelle nodes and the component monomers converged to the same value at long times (t ≈ 300). It is usual to normalize the measured D values with a limiting zero density value, D0. The concentration-dependent aggregation number found in the simulations presents a problem for the determination of Dm/D0 because D0 will depend on m. The diffusion coefficient of the whole polymer was calculated. (Recall that the polymer is represented by a simple blob rather than a necklace of beads as is more usual, so the procedure was therefore straightforward.) The mean-square displacement was that of the center of mass of the polymer. When micelle aggregation number is constant, then D0 can be determined by extapolation to φ f 0. With varying aggregation number, however, the hydrodynamic radius of the micelles, rhm, will be different at each packing fraction, and hence D0 is φ-dependent. Employing the simple approximation that the micelle is like a star polymer with m arms, the radius of gyration of the micelle, rgyr, is approximately (3 - (2/m))ν in reduced units, where ν is the Flory exponent.52 For a theta solvent, ν ) 1/2, and for a good solvent, ν ) 3/5. Furthermore, the ratio of the hydrodynamic radius to the radius of gyration (rhm/rgyr) increases with m,52 which can also depend on solvent quality. Using the Stokes-Einstein relationship, this gives, for example, D0 ) 4.36 × 10-2 for m ) 4 and D0 ) 3.69 × 10-2 for m ) 12 in a θ solvent. Taking D0 to be 1/(6π x3-(2/m)) (i.e., assuming that rhm = rgyr ) (3 - (2/m))3/5), D/D0 is plotted as a function of φ in Figure 12. In the Figure, D/D0 shows a near-logarithmic dependence on packing fraction in the range considered, which is consistent with experiment at intermediate concentrations.25 At very low concentrations, a sudden increase in the measured aggregation number (anomalous micellization) is often observed,25 the cause of which is uncertain. This feature is not evident in these simulations.
Cass et al.
Figure 13. Simulation viscosity divided by packing fraction plotted against packing fraction. By analogy with the standard Huggins plot of experimental data (i.e., ηp/c vs c), the data reduction is convenient for determining the Huggins coefficient by linear regression (dashed line) to the data points that are shown as circles. N ) 1000, Ea ) 16, and ra ) 6.0. φ ) 1.85φ′.
As discussed in section 2, the present level of coarse graining eliminates most of the degrees of freedom responsible for the intrinsic viscosity, [η]. In fact, the intrinsic viscosity is completely absent for the current model for nonassociating polymers, which cannot form nodes, because these states evolve solely with effective intermolecular excluded volume interactions as given by the potential of eq 1. When polymers associate with nodes, the situation becomes more complex because the node-polymer bead interaction represents some of the internal degrees of freedom of the polymer, hence part of the intrinsic viscosity is then an outcome of the simulation. The node-core interaction is included in the stress, the stress correlation function in eq 6, and hence the viscosity. There is a partial intrinsic viscosity contribution to the computed viscosity, and the treatment of the data should now use ηp rather than ηip. Simulation gives
ηp ) [η]′c + kH[η]2c2 + ...
(30)
where [η]′ is part of the full intrinsic viscosity [η], which can be written as [η]′ ) f [η] in which f is a dimensionless quantity between 0 and 1. When the relation φ ) (2/5)[η]c is used, eq 30 can be rewritten as
ηp ) 6.25kHφ + 2.5f φ
(31)
The f term should increase with the fraction of polymer molecules associated with nodes at any given time (for nonassociating polymers f ) 0). Given that even at low concentration polymers are almost always associated with a node, it is reasonable to assume that f is constant to first approximation. Figure 13 shows the viscosity found by simulation over φ as a function of φ and a linear fit to this data. The fit gives f ) 0.75 and kH ) 5.75. Although there is only one new degree of freedom on the mesoscale, the node-core interaction represents many degrees of freedom on the atomistic and even Kuhn length scales, so the relatively large magnitude of f is perhaps not surprising. Quintana and co-workers found kH ≈ 0.75 for nearly symmetric diblock polymers in selective solvents.56 This would be slightly in excess of the value for a homopolymer in a theta solvent but not as high (56) Quintana, J. R.; Villacampa, M.; Katime, I. A. Macromolecules 1993, 26, 606.
Brownian Dynamics Simulations of Diblock Copolymers
in magnitude as the value found in the present study or the experimental values for telechelics with large n-alkyl end caps.57,58 The fact that their data were described well by the Huggins equation was interpreted in terms of the micellar dimensions being invariant over the concentration range studied. The value of f will in fact not be constant but will depend on the fraction of polymers associated (and hence experiencing a node-polymer interaction) and the balance of intra- and interpolymer contributions to the node-polymer potential (which will depend on the separation of polymers from their nodes). These factors will lead to a more variable f at lower packing fractions; in general, as φ tends to zero, f tends to zero because the fraction of polymers associated with nodes diminishes.
IV. Conclusions The main novelty of this work is that we have invented a new computer model for associative polymers that accounts for the two widely separated time scales dominating the physical behavior of diblock polymers. The Born-Oppenheimer approximation to relax the electronic degrees of freedom of atoms and molecules essentially instantaneously on the time scale of the relaxation of the nuclei has been widely applied in condensed matter simulations to include explicit quantum effects “on the fly” (e.g., see ref 59). The proposed coarse-grained computer simulation model of associative polymers has much in common with this approximation in that we assume the micelle core instantaneously follows the positions of the polymer cores (i.e., on their time scale). This approximation is valid because it reflects the wide separation of the two time scales (core center of mass and hydrophobe disengagement/re-engagement) present in the real systems. Here we demonstrate the proof of concept and that the trends modeled agree with experiment (e.g., see Figure 10 for the concentration dependence of the aggregation numbers). In future studies, it would be possible to improve agreement with experiment for particular chemical systems by fine tuning the parameters. The micellization of a coarse-grained soft blob model for diblock coblocks in solution has been investigated using Brownian dynamics simulations. The current model applies to the limit of a large hydrophilic group (represented by a single soft-core particle with a radius on the order of the radius of gyration of the polymer) attached to a much smaller hydrophobic block. It is best applied in the dilute and the initial semidilute regimes prior to entanglement effects. This coarse graining allows a large number of polymers to be simulated over dynamically and rheologically important time scales that cover polymer association/dissociation kinetics from micelles and significant diffusion and interdiffusion (57) Jenkins, R. D. Ph.D. Thesis, Lehigh University, 1990. (58) Hough, R.-L.; Ph.D. Thesis, University of Wales, 1993. (59) Carr, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471.
Langmuir, Vol. 23, No. 12, 2007 6587
of the micelles. Despite its simplicity, the results of the proposed model compare favorably with experimental trends and help identify the origin of experimental behavior in terms of key length scales and interaction energies on the mesoscale level of description used here. The micelle size distribution is similar to that found in a previous simulation study of associating coblock polymers in solution.54 The aggregation numbers found are of the same order as those found experimentally for long-chain polymers with small hydrophobes.25 The concentration at which the micelles start to form, the critical micelle concentration, is typical of coblock polymers, being very low compared to that of typical surfactant systems. Our simulation systems formed micelles at all of the concentrations explored. The dependence of aggregation number on packing fraction fits well with a simplified version of Leibler et al.’s theory.47 This dependence is weak for the case of large aggregation numbers. However, with the small micelles, which result from highly asymmetric block copolymers, there is a significant impact on aggregation number over the same packing fraction range. The packing fraction dependence of the selfdiffusion coefficient at intermediate concentrations mirrors that found in experiment and exhibits a decrease as a function of the logarithm of the packing fraction. The Huggins coefficient can only be approximated (kH ≈ 6) because of the incomplete inclusion of intrinsic viscosity in the viscosity used in the simulation model. The simulations well reproduce the low-to-intermediate packing fraction properties of solutions of associating asymmetric block copolymers but probably exceed their bounds of validity for much higher concentrations than the overlap concentration (for packing values of ∼0.5 here). The model should be readily extendible to the case of associating polymers with multiple stickers for which more extensive intrapolymer interactions and micellar structures are possible. The model could also be used in other coarse-grained simulation methods in which the effective interactions are typically rather soft, such as dissipative particle dynamics (DPD).60,61 Unlike real colloidal liquids, the micelles are not static structures because the free diblock molecules in solution are in dynamic equilibrium with those in the micelle, so it would be less well founded to model the whole micelle as a single blob. Acknowledgment. We thank the Engineering and Physical Sciences Research Council of Great Britain (EPSRC) for funding this project. LA063210J (60) Hoogerbrugge, P. J.; Koelmann, J. M. V. A. Europhys. Lett. 1992, 19, 155. (61) Nol, P. E.; Warren, P. B. Europhys. Lett. 1995, 30 191.
