4076
Langmuir 2006, 22, 4076-4083
Monte Carlo Simulations of Micellization in Model Ionic Surfactants: Application to Sodium Dodecyl Sulfate Daniel W. Cheong and Athanassios Z. Panagiotopoulos* Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544 ReceiVed December 29, 2005. In Final Form: February 2, 2006 A lattice model for ionic surfactants with explicit counterions is proposed for which the micellization behavior can be accurately determined from grand canonical Monte Carlo simulations. The model is characterized by a few parameters that can be adjusted to represent various linear surfactants with ionic headgroups. The model parameters have a clear physical interpretation and can be obtained from experimental data unrelated to micellization, namely, geometric information and solubilities of tail segments. As a specific example, parameter values for sodium dodecyl sulfate were obtained by optimizing for the solubility of hydrocarbons in water and the structural properties of dodecane. The critical micelle concentration (cmc), average aggregation number, degree of counterion binding, and their dependence on temperature were determined from histogram reweighting grand canonical Monte Carlo simulations and were compared to experimental results. The model gives the correct trend and order of magnitude for all quantities but underpredicts the cmc and aggregation number. We suggest ways to modify the model that may improve agreement with experimental values.
1. Introduction It is widely known that surfactants with the right balance of hydrophobic and hydrophilic moieties self-assemble into micelles in solution. The concentration at which micellization occurs, known as the critical micelle concentration (cmc), depends largely on the surfactant architecture, temperature, type of solvent, and the presence of other solutes such as salt or cosurfactants.1-4 Self-assembly and micellization of surfactants are important in numerous practical applications. For instance, surfactants are important for their cleaning properties and are the main ingredients in many common household cleaning products.2,4 Surfactant selfassembly is also important in nanotechnology, where materials can be induced to organize spontaneously into suitable templates to form nanostructures.5,6 The principle of self-assembly is also relevant in understanding the behavior of more complex structures such as biomembranes.1 Surfactant micellization is a relatively well-studied topic. The structure and micellization properties of various surfactants have been investigated through experimental techniques such as light scattering, fluorescence quenching, calorimetry, potentiometry, and so forth.7-17 Values of the critical micelle concentration and * Corresponding author. E-mail:
[email protected]. (1) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: Orlando, FL, 1985. (2) Rosen, M. J. Surfactants and Interfacial Phenomena, 2nd ed.; Wiley: New York, 1989. (3) Tanford, C. The Hydrophobic Effect: Formation of Micelles and Biological Membranes, 2nd ed.; Krieger Publishing Co.: Malabar, FL, 1991. (4) Myers, D. Surfactant Science and Technology, 2nd ed.; VCH Publishers: New York, 1992. (5) McGrath, K. M.; Dabbs, D. M.; Yao, N.; Aksay, I. A.; Gruner, S. M. Science 1997, 277, 552. (6) Faul, C. F. J.; Antonietti, M. AdV. Mater. 2003, 15, 673. (7) Hayashi, S.; Ikeda, S. J. Phys. Chem. 1980, 84, 744. (8) Bezzobotnov, V. Y.; Borbely, S.; Cser, L.; Farago´, B.; Gladkih, I. A.; Ostanevich, Y. M.; Vass, S. J. Phys. Chem. 1988, 92, 5738. (9) Okazaki, M.; Polyakov, N. E.; Toriyama, K. J. Chem. Phys. 1995, 99, 6452. (10) Paula, S.; Sues, W.; Tuchtenhagen, J.; Blume, A. J. Phys. Chem. 1995, 99, 11742. (11) Bijma, K.; Engberts, J. B. F. N. Langmuir 1997, 13, 4843. (12) Bales, B. L.; Messina, L.; Vidal, A.; Peric, M.; Nascimento, O. R. J. Chem. Phys. 1998, 102, 10347. (13) McMahon, C. A.; Hawrylak, B.; Marangoni, D. G.; Palepu, R. Langmuir 1999, 15, 429.
average aggregation number for a wide range of surfactants can be readily found in the literature.2 Theories have also been developed to predict the micellization properties of surfactants.18-27 More recently, computer simulations have been increasingly used to study the structure and thermodynamics of micelles.28-60 (14) Kira´ly, Z.; Deka´ny, I. J. Colloid Interface Sci. 2001, 242, 214. (15) Ranganathan, R.; Peric, M.; Medina, R.; Garcia, U.; Bales, B. Langmuir 2001, 17, 6765. (16) Chatterjee, A.; Moulik, S. P.; Sanyal, S. K.; Mishra, B. K.; Puri, P. M. J. Phys. Chem. B 2001, 105, 12823. (17) Aswal, V. K.; Goyal, P. S.; Amenitsch, H.; Bernstorff, S. Pramana 2004, 63, 333. (18) Ruckenstein, E.; Nagarajan, R. J. Phys. Chem. 1981, 85, 3010. (19) Dill, K. A.; Koppel, D. E.; Cantor, R. S.; Dill, J. D. Nature 1984, 309, 42. (20) Ben-Shaul, A.; Szleifer. I.; Gelbart, W. M. J. Chem. Phys. 1985, 83, 3597. (21) Szleifer, I.; Ben-Shaul, A.; Gelbart, W. M. J. Chem. Phys. 1986, 85, 5345. (22) Nagarajan, R.; Wang, C. C. Langmuir 1995, 11, 4673. (23) Zoeller, N.; Lue, L.; Blankschtein, D. Langmuir 1997, 13, 5258. (24) Guerin, C. B. E.; Szleifer, I. Langmuir 1999, 15, 7901. (25) Moroi, Y.; Nobuyuki, Y. Langmuir 1997, 13, 3909. (26) Srinivasan, V.; Blankschtein, D. Langmuir 2003, 19, 9932. (27) Srinivasan, V.; Blankschtein, D. Langmuir 2003, 19, 9946. (28) Jo¨nsson, B.; Edholm, O.; Teleman, O. J. Chem. Phys. 1986, 85, 2259. (29) Watanabe, K.; Ferrario, M.; Klein, M. L. J. Chem. Phys. 1988, 92, 819. (30) Karaborni, S.; O’Connell, J. P. Langmuir 1990, 6, 905. (31) Karaborni, S.; O’Connell, J. P. J. Phys. Chem. 1990, 94, 2624. (32) Maillet, J. B.; Lachet, V.; Coveny, P. V. Phys. Chem. Chem. Phys. 1999, 1, 5277. (33) Tieleman, D. P.; van der Spoel, D.; Berendsen, H. J. C. J. Phys. Chem. B 2000, 104, 6380. (34) Marrink, S. J.; Tieleman, D. P.; Mark, A. E. J. Phys. Chem. B 2000, 104, 12165. (35) Bogusz, S.; Venable, R. M.; Pastor, R. W. J. Phys. Chem. B 2000, 104, 5462. (36) Bogusz, S.; Venable, R. M.; Pastor, R. W. J. Phys. Chem. B 2001, 105, 8312. (37) Garde, S.; Yang, L.; Dordick, J. S.; Paulaitis, M. E. Mol. Phys. 2002, 100, 2299. (38) Shelley, J.; Watanabe, K.; Klein, M. L. Int. J. Quantum Chem. 1990, 38, 103. (39) MacKerell, A. D. J. Phys. Chem. 1995, 99, 1846. (40) Schweighofer, K. J.; Essmann, U.; Berkowitz, M. J. Phys. Chem. B 1997, 101, 3793. (41) Bruce, C. D.; Berkowitz, M. L.; Perera, L.; Forbes, M. D. E. J. Phys. Chem. B 2002, 106, 3788. (42) Bruce, C. D.; Senapati, S.; Berkowitz, M. L.; Perera, L.; Forbes, M. D. E. J. Phys. Chem. B 2002, 106, 10902. (43) Karaborni, S.; van Os, N.; Esselink, K.; Hilbers, P. Langmuir 1993, 9, 1175. (44) Smit, B.; Esselink, K.; Hilbers, P.; van Os, N. Langmuir 1993, 9, 9. (45) Wijmans, C. M.; Linse, P. Langmuir 1995, 11, 3748.
10.1021/la053511d CCC: $33.50 © 2006 American Chemical Society Published on Web 03/09/2006
Micellization in Model Ionic Surfactants
There are two main classes of simulation models that can be used to study surfactant systems: detailed atomistic models and coarse-grained models.61,62 A detailed atomistic model consists of a large number of water and surfactant molecules, and all water-water, water-surfactant, and surfactant-surfactant interactions are considered explicitly. All possible separations, orientations, and conformations are also considered. In theory, this atomistic approach would be the most rigorous and accurate approach. However, the characteristic time for micelles to form in solution and come to equilibrium is generally in the microsecond range and above61,62 and thus is much larger than the time scales currently accessible by atomistic simulations.61,62 This makes atomistic simulations unfeasible for studying spontaneous micelle formation and equilibrium properties. However this approach is useful in determining more detailed and accurate information about the local structure of surfactant solutions. The size and shape, short-time dynamics, and counterion condensation of detailed atomistic models of various surfactant micelles have been investigated by molecular dynamics simulations.28-42 An alternative to atomistically detailed models is coarse-grained models in which a number of atoms are grouped together. These models are usually able to describe micellization although they lack the detailed structural information that the atomistic model can provide. These models can in turn be divided into continuum and lattice models. Continuum models are more realistic, but lattice models are more computationally efficient. Numerous simulation studies involving coarse-grained continuum43-49 and lattice50-60 models have been performed. Lattice models have been used extensively to study nonionic micelles using Monte Carlo simulations.50-60 The most common lattice model of nonionic surfactants is probably that of Larson.51 This model places the surfactant molecules on a fully occupied cubic lattice and restricts the interactions as well as the direction of growth to the 26 neighboring sites. Care52 used a very similar lattice model but considers only adjacent sites (no diagonal sites) for chain growth and interactions. Using these lattice models, Monte Carlo simulations were performed, and the thermodynamic properties50-56 and the size and shape distributions58-60 of the micelles have been determined. Ionic surfactants have been studied much less than their nonionic counterparts. Bhattacharya et al.63,64 used Monte Carlo simulations to study a 2D continuum model of ionic surfactants. No explicit counterions were considered. Interaction potentials (46) Palmer, B. J.; Liu, J. Langmuir 1996, 12, 746. (47) von Gottberg, F.; Smith, K.; Hatton, T. J. Chem. Phys. 1997, 106, 9850. (48) Viduna, D.; Milchev, A.; Binder, K. Macromol. Theory Simul. 1998, 7, 649. (49) Milchev, A.; Bhattacharya, A.; Binder, K. Macromolecules 2001, 34, 1881. (50) Haan, S.; Pratt, L. Chem. Phys. Lett. 1981, 79, 436. (51) Larson, R.; Scriven, L.; Davis, H. J. Chem. Phys. 1985, 83, 2411. (52) Care, C. J. Phys. C: Solid State Phys. 1987, 20, 689. (53) Brindle, D.; Care, C. J. Chem. Soc., Faraday Trans. 1992, 88, 2163. (54) Bernardes, A.; Henriques, V.; Bisch, P. J. Chem. Phys. 1994, 101, 645. (55) Floriano, M. A.; Caponetti, E.; Panagiotopoulos, A. Z. Langmuir 1999, 15, 3143. (56) Panagiotopoulos, A. Z.; Floriano, M. A.; Kumar, S. K. Langmuir 2002, 18, 2940. (57) Care, C. J. Chem. Soc., Faraday Trans. 1987, 83, 2905. (58) Larson, R. J. Chem. Phys. 1992, 96, 7904. (59) Mackie, A. D.; Panagiotopoulos, A. Z.; Szleifer, I. Langmuir 1997, 13, 5022. (60) Nelson, P. H.; Rutledge, G. C.; Hatton, T. A. J. Chem. Phys. 1997, 107, 10777. (61) Gelbart, W. M.; BenShaul, A. J. Phys. Chem. 1996, 100, 13169. (62) Shelley, J. C.; Shelley, M. Y. Curr. Opin. Colloid Interface Sci. 2000, 5, 101. (63) Bhattacharya, A.; Mahanti, S.; Chakrabarti, A. J. Chem. Phys. 1998, 108, 10281. (64) Bhattacharya, A.; Mahanti, S. J. Phys.: Condens. Matter 2001, 13, 1413.
Langmuir, Vol. 22, No. 9, 2006 4077
were described via a Lennard-Jones potential, a bond-bending potential, and a screened Coulombic potential. Spontaneous micellization was observed in the simulations, but no quantitative predictions of bulk micellization properties could be made because the model was 2D. Other simulation studies involving ionic surfactants focus on the interaction between like-charged colloidal particles (i.e., charged micelles)65,66 or the interaction and complexation of polyelectrolytes with charged micelles.67-69 In such studies, micelles of a certain size and charge are assumed to exist in the system already, and the simulations attempt to study their behavior in solution. An overview of the use of computer simulations to study surfactants can be found in several review articles.62,70,71 To our knowledge, no systematic simulation study of the micellization of ionic surfactants in three dimensions has been performed. Simulations have also not been used to predict micellization properties, such as the critical micelle concentration, size and shape of micelles, degree of counterion condensation, and so forth. In this work, we have developed a simple lattice model that can be used to obtain quantitative predictions of the micellization properties of ionic surfactants. As an example, we studied sodium dodecyl sulfate (SDS) and compared the predicted properties with experimental results to evaluate the accuracy of the model. The model can also be used to study other linear surfactants by adjusting its parameters, which have a clear physical significance and can be determined from experimental data unrelated to micellization (e.g., bulk solubilities). The rest of this article is structured as follows. In section 2, we describe the model that we have developed and used to represent the ionic surfactant and describe how parameters were obtained to model SDS. In section 3, we describe some of the simulation techniques used in this work. We describe how we obtain the cmc from our simulations, and the results for the cmc are presented and compared with experimental values in section 4. The micelle size and degree of counterion binding are presented in section 5. Finally, we conclude by summarizing our results and discussing the limitations of our current model in section 6. We also suggest possible improvements to the model to obtain better quantitative agreement with experiment.
2. Model In our model, a surfactant molecule is represented as a chain of sites on a simple cubic lattice consisting of Nh headgroups of diameter σh and Nt tail groups of diameter σ. Each headgroup carries a charge of Qh and is neutralized by a number of monovalent counterions to balance the charge of the headgroup exactly. All counterions have the same diameter σ as a tail group, which also defines the lattice spacing. The headgroup diameter σh can be made arbitrarily large relative to σ. Tail group connectivity is limited to the nearest neighbors, resulting in a coordination number ztt ) 6 in three dimensions (Care model52). This makes the chain more rigid, and bond angles are restricted to either 90 or 180°. This choice is justified later in this section. The connectivity between a headgroup and a tail group is limited to the number of lattice sites that satisfy the constraint
σht e dht < (σht + 1)
(1)
(65) Linse, P.; Lobaskin, V. J. Chem. Phys. 2000, 112, 3917. (66) Linse, P. Philos. Trans. R. Soc. London, Ser. A 2001, 359, 853. (67) Linse, P.; Wallin, T. J. Phys. Chem. B 1997, 101, 5506. (68) Wallin, T.; Linse, P. Langmuir 1998, 14, 2940. (69) Hansson, P. Langmuir 2001, 17, 4167. (70) Karaborni, S.; Smit, B. Curr. Opin. Colloid Interface Sci. 1996, 1, 411. (71) Larson, R. Curr. Opin. Colloid Interface Sci. 1997, 2, 361.
4078 Langmuir, Vol. 22, No. 9, 2006
Cheong and Panagiotopoulos
Figure 1. Schematic diagram of ionic surfactants and their counterions. The diagram is a 2D illustration of the 3D molecules simulated in this work. Table 1. Parameters Used to Model Sodium Dodecyl Sulfate parameter
value
Nh Nt Qh σ σh tt/kB
1 6 -1 3.97 Å 5.96 Å -600 K
where σht ) (σ + σh)/2 is the collision diameter between a headgroup and a tail group and dht is the distance between a headgroup and a tail group. Similarly, if there are multiple headgroups, then connectivity between two headgroups is limited to the sites that are within one lattice spacing of the headgroup diameter. All counterions are free to be on any unoccupied position on the lattice. A schematic of the model is given in Figure 1. The charged groups interact via the Coulombic potential and a hard core repulsion. Headgroups and counterions thus interact via
{
qiqje2 1 for rij g (σi + σj) 2 4πDD r 0 ij Uij(rij) ) 1 for rij < (σi + σj) +∞ 2
(2)
where qi and qj are the charges of the charged particles, e is the elementary charge, D and D0 are the dielectric constant of the intervening medium and the dielectric permittivity of vacuum, respectively, and rij is the distance between the two charged particles. Interactions between tail groups are of strength tt and are limited to the same range as the tail-tail connectivity of six nearest neighbors. No tail-tail interactions occur at greater distances. A hard core repulsion also exists between a tail group and a headgroup according to the second condition of eq 2. Only a single tail group or counterion can occupy a lattice site, enforcing an excluded volume condition for the tail-counterion and tailtail interactions. No short-range interaction besides volume exclusion exists between headgroups or between the head and tail groups. The parameters in the model can be chosen to model different surfactants. To model SDS, the following parameters were used and are listed in Table 1. The reasons behind these choices for the parameters are explained in the following paragraphs. The lattice spacing was set to the value σ ) 3.97 Å. A tail group was set to represent 2 methyl groups, resulting
Figure 2. Solubility of hydrocarbons in water as a function of temperature. The solid lines are the low-density branches of the phase diagrams, and the filled circles are experimental values.73,74 Red corresponds to T3 and C6, blue to T4 and C8, orange to T5 and C10, and pink to T6 and C12. With appropriate parameters, the solubilities of hydrocarbons in water at 300 K are matched rather well.
in Nt ) 6 tail groups for the 12 methyl groups in SDS. A single headgroup (Nh ) 1) of diameter σh ) 1.5σ ) 5.96 Å and charge Qh ) -1 was used to represent the sulfate group. A single lattice site of charge +1 was used to represent the sodium Na+ counterion. The tail-tail interaction strength was set to be tt/kB ) -600 K, resulting in attractive interactions for nearest-neighbor tail-tail contacts. A dielectric constant of D ) 80 was chosen to model an aqueous surfactant solution around room temperature. On first observation, it may seem natural to have each tail group represent a single methyl group and set the diameter to be the van der Waals diameter of the methyl group σvdw ≈ 3.4 Å. However, because the tail beads are modeled as single lattice sites, this will overestimate the end-to-end distance of the hydrocarbon chain where the average C-C bond length is only 1.54 Å. Furthermore, because we have not included any torsional potential in our model, it will not account for the rigidity present in real hydrocarbon chains. Allowing a tail group to represent more than one methyl group seems to be a reasonable compromise to account for the hydrocarbon chain rigidity and the length of the C-C bond. We also decided to use a coordination number of z ) 6 to define the chain connectivity. Allowing only nearestneighbor connectivity reduces the end-to-end distance and increases the rigidity of the chain. This is in contrast to previous studies of nonionic surfactants in our group,55,56 where z ) 26 was used. We obtained the phase envelopes for homopolymer chains of Nt ) 3, 4, 5, and 6 tail beads, similar to an earlier study of longer homopolymer chains by Panagiotopoulos et al.72 The low-density branch of the phase envelope corresponds to the solubility of the chains. Figure 2 shows the reported solubility of various hydrocarbons in water from the literature73,74 together with the low-density branch of the phase diagram for T3-T6, where Tx denotes a chain with x tail beads. The simulations cannot predict the nonmonotonic temperature dependence of the solubilities because of the lack of explicit solvent. However, because we are interested in studying the spontaneous micellization of surfactants at room temperature, we decided to match (72) Panagiotopoulos, A. Z.; Wong, V.; Floriano, M. A. Macromolecules 1998, 31, 912. (73) Ma¸ czyn´ski, B.; Wis´niewska-Gocłowska, B.; Go´ral, M. J. Phys. Chem. Ref. Data 2004, 33, 549. (74) Sutton, C.; Calder, J. A. EnViron. Sci. Technol. 1974, 8, 654.
Micellization in Model Ionic Surfactants
the solubility data at 300 K. Figure 2 demonstrates that by allowing each tail bead to correspond to two methyl groups, each empty lattice site to correspond to one water molecule, and by using the appropriate value for the interaction strength tt the solubility of hydrocarbons in water for the range of C6-C12 at T ) 300 K can be matched reasonably well. The values to obtain this matching, Nt ) 6 and tt/kB ) -600 K, were adopted in this work. There is a narrow range of values of tt/kB for which the solubilities can be reasonably matched, but the effect on the calculated micellization properties should be minimal. The size of a tail bead, which also defines the lattice spacing, must be determined as well. This was done by matching the liquid density of C12 at 298 K, which can be obtained from standard handbooks75 to be F ) 0.751 kg/m3. By matching the liquid coexistence density of T6 at 298 K to the experimental liquid density of dodecane, a diameter of σ ) 3.97 Å was obtained. To verify that this value of σ will reasonably depict the structure of dodecane, the end-to-end distance 〈r2〉1/2 and radius of gyration 〈s2〉1/2 were compared to the values from previous atomistic simulations.76 With σ ) 3.97 Å, our simulations yielded 〈r2〉1/2 ) 10.17 Å and 〈s2〉1/2 ) 4.35 Å. Although slightly overestimated when compared to 〈r2〉1/2 ) 9.79 Å and 〈s2〉1/2 ) 3.52 Å, obtained from atomistic simulations,76 the structural properties of dodecane are reasonably reproduced. The end-to-end distance and radius of gyration were also calculated for chains with z ) 26 and were found to be even higher than the values at z ) 6. This reinforced the choice of z ) 6. The bond length and the van der Waals radius of the methyl group cannot be simultaneously matched using a hard core potential, and the bond angles in a real hydrocarbon chain also cannot be reproduced on a lattice. Consequently, it is not possible to match the hydrocarbon structure exactly using our simple lattice model. The size of the sulfate headgroup was determined by accounting for the sulfur-oxygen bond length and the van der Waals radii for the oxygen and sulfur. With an S-O bond length of 1.57 Å and a van der Waals radius for the oxygen of 1.4 Å, we have approximated the diameter of the sulfate group to be approximately twice the length of the S-O bond and oxygen van der Waals radius or 1.5 times the lattice spacing, σh ≈ 5.96 Å ) 1.5σ. Only the Nh ) 1 headgroup is necessary for the sulfate group, which carries a charge of Qh ) -1 corresponding to the charge of the sulfate group. For simplicity, the counterions also occupy a single lattice site, giving them a diameter of 3.97 Å, which is larger than the ionic diameter that is around 2 Å. Clearly, the lattice model cannot reproduce all of the surfactant’s properties. Because the hydrophobic interactions are the main driving force of micellization, we decided that it was most important to describe the behavior of the hydrocarbon tail as well as possible. Although the chain structure is greatly simplified, this set of parameters provides a balance to reproduce as closely as possible the solubility of hydrocarbons in water, the liquid density of dodecane, and the end-to-end distance and radius of gyration of dodecane.
3. Simulation Method Grand canonical Monte Carlo (gcmc) simulations were performed in cubic boxes of length L under periodic boundary conditions. A typical mix of 50% insertion/deletion moves, 35% single chain displacements, 5% cluster displacements, and 10% single counterion displacements was performed. To handle the long-ranged Coulombic potential, the Ewald summation method (75) Perry, R. H.; Green, D. W., Eds. Perry’s Chemical Engineering Handbook, 6th ed.; McGraw-Hill: New York, 1984. (76) Avitabile, G.; Tuzi, A. J. Polym. Sci.: Polym. Phys. Ed. 1983, 21, 2379.
Langmuir, Vol. 22, No. 9, 2006 4079
was used with conducting boundary conditions, 518 Fourierspace wave vectors, and real-space damping parameter κ ) 5. Because we are performing lattice simulations, all of the Coulombic interactions can be precomputed and stored in an array. Subsequent energy calculations are reduced to a table lookup procedure, speeding up the simulations significantly.77 The chain insertions and deletions were performed using a configurational-bias method78-80 that allows chain growth only into unoccupied spaces as described in previous studies of homopolymer chains72 and nonionic surfactants.55 During an insertion, a headgroup is first inserted into an available lattice site, and a counterion is subsequently inserted into an unoccupied lattice site with a distance bias.81 The chain of tail beads is then grown beginning from an available neighboring site of the inserted headgroup. In the single chain and counterion displacement move, a chain or counterion is selected randomly and displaced in a random direction. Cluster displacements were performed according to the convention described by Mackie et al.59 Two surfactant molecules are considered to be a part of the same cluster if a tail group from one surfactant is in one of the 26 nearest- and next-nearest-neighbor sites of a tail group from the other surfactant. Starting from a randomly selected chain, a cluster is defined and displaced in a random direction. The move is accepted according to the normal Metropolis acceptance criteria. However, if a cluster is moved such that the number of chains in the cluster changes, that move has to be rejected to satisfy detailed balance. Using these methods, simulations were performed in boxes of length L ) 79.4, 119.1, and 158.8 Å at varying temperatures. Three sets of independent runs using different random number sequences were performed under the same thermodynamic conditions to obtain statistical uncertainties. Histograms from the different runs are combined using multihistogram reweighting techniques.82 As described in the work on nonionic surfactants,55 this technique allows the prediction of the cmc at varying temperatures to be obtained quite efficiently. The results for the cmc’s are presented and compared with experimental values in section 4. During the simulations, we also keep track of the micelle cluster size and the counterion distribution around a micelle. These results are presented in section 5.
4. Determination of the Critical Micelle Concentration The cmc is a key surfactant property. It represents the concentration of surfactants in solution necessary for micellization to occur. Micellization can usually be recognized by a sharp break in the solution bulk properties at the cmc. In this work, we have obtained the cmc from the osmotic pressure curve according to the method described by Floriano et al.55 This procedure is described below. The osmotic pressure can be determined from the simulations by the relationship ln Ξ ) βPV, where Ξ is the grand partition function. We performed a series of gcmc simulations over a range of chemical potentials and temperatures, ensuring sufficient overlaps of the surfactant densities between each successive run. The results were then combined through histogram reweighting. At low chemical potentials, the surfactant molecules exist as monomers. Upon increasing the chemical potential, the surfactant (77) Panagiotopoulos, A. Z.; Kumar, S. K. Phys. ReV. Lett. 1999, 83, 2981. (78) Frenkel, D.; Mooij, G. C. A. M.; Smit, B. J. Phys.: Condens. Matter 1992, 4, 3053. (79) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59. (80) de Pablo, J. J.; Laso, M.; Siepmann, J. I.; Suter, U. W. Mol. Phys. 1993, 80, 55. (81) Orkoulas, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 1994, 101, 1452. (82) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1988, 61, 2635.
4080 Langmuir, Vol. 22, No. 9, 2006
Cheong and Panagiotopoulos Table 2. Critical Micelle Concentration (in mM) for Model SDS from Simulations in Boxes of L ) 79.4, 119.1, and 158.8 Åa T(K)
L ) 79.4 Å
L ) 119.1 Å
L ) 158.8 Å
400 375 350 325 300
10.6(3) 7.3(2) 4.79(2) 2.78(2) 1.21(2)
9.22(9) 5.79(1) 3.47(2) 2.09(2) 1.06(3)
8.3(2) 5.4(1) 3.10(8) 1.64(5) 0.97(1)
a Statistical uncertainties in parentheses refer to the last decimal place shown.
Figure 3. Snapshot of two SDS micelles in a simulation box of L ) 158.8 Å at T ) 300 K. Blue beads are the headgroups, and gray beads are the tail groups. Red beads are counterions that are condensed, and green beads are free counterions. The micelles have an aggregation number of M ≈ 40 at an overall concentration of about 36 mM.
Figure 4. Simulation results of the osmotic pressure P(kPa) vs the surfactant concentration (mM) at different temperatures for L ) 119.1 Å. From top to bottom, the curves are for T ) 375, 350, 325, and 300 K. The breaks in the curves correspond to the onset of micellization. The cmc is defined as the concentration at the inflection point of the curve.
molecules begin to aggregate and form micelles. The micelles are fairly spherical with a well-defined and densely packed core of tail groups surrounded by charged headgroups. A snapshot of a system with two micelles is shown in Figure 3. To determine the cmc, we determined the osmotic pressure P and plotted it against the surfactant concentration at varying temperatures as shown in Figure 4. We observe that there are two distinct regions of different slopes in the curves. At low concentrations, the osmotic pressure curve follows a higher slope that corresponds to the single surfactant molecules in solution. At some concentration that depends on temperature, there is a break in the curve, and a second straight line with a lower slope continues at higher concentrations. Thermodynamically, this indicates that above a certain concentration surfactant molecules in the system begin forming aggregates and the concentration of free surfactants will remain roughly constant around the cmc.
Figure 5. Critical micelle concentration of SDS as a function of temperature. Open symbols are experimental values: diamonds are from Blume et al.;10 squares are from Goddard and Benson;84 and circles and triangles are from Chatterjee et al.16 Filled circles are the simulation results for L ) 158.8 Å. Error bars are smaller than the symbol size. Lines are guides for the eye.
We define the cmc as the concentration at the point of inflection of the osmotic pressure curve. This corresponds to the concentration at which micelles begin to appear in the system. At higher temperatures, the break becomes increasingly gradual, and the cmc becomes harder to define accurately. This effect is completely analogous to the difficulties of determining the cmc experimentally at higher temperatures.83 Following this definition, we obtained the cmc of SDS at L ) 79.4, 119.1, and 158.8 Å. The results of the cmc (mM) as a function of temperature are given in Table 2. The cmc curve for the largest box size considered (L ) 158.8 Å) is shown in Figure 5 together with the experimental results obtained by calorimetric, conductometric, and potentiometric methods.10,16,84 The model correctly predicts the order of magnitude of the cmc and the general trend where, above room temperature, the cmc increases with an increase in temperature. However, the simulation results consistently predict too low a cmc when compared to the experimental results. Because micellization is primarily driven by the hydrophobic effect, the low cmc suggests that the tailtail attractive interactions are too strong or that there are too many favorable interactions occurring. The value of the interaction strength tt/kB could be adjusted to address this issue. The effect of the headgroup (hydrophilic interactions) is usually small in comparison to the hydrophobic interactions of the hydrocarbon tail. However, increasing the charge of the headgroup will increase the electrostatic repulsion between them, inhibiting micellization and raising the cmc. A charge of Qh ) -1 seems to be the most appropriate for SDS, but our model does not include short-range repulsions between headgroups and between the tail and headgroups. Short-range repulsions between headgroups will (83) Emerson, M. F.; Holtzer, A. J. Phys. Chem. 1967, 71, 3320. (84) Goddard, E.; Benson, G. Can. J. Chem. 1957, 35, 986.
Micellization in Model Ionic Surfactants
also lead to an increase in the cmc for the same reasons as increasing the headgroup charge. Therefore, the inclusion of a short-range repulsion between headgroups could play an important role and should be included in future studies using this model. The experimental cmc curves exhibit a minimum around 300 K. At low temperatures, hydrophobic interactions between the surfactant tail and the water molecules cause distortions in the local water structure around the surfactant. Micellization is strongly favored because it minimizes this distortion in the water structure. At the same time, electrostatic repulsions between the like-charged headgroups work against micellization, which will bring them closer to each other. As the temperature is raised, the electrostatic repulsions between the headgroups become weaker, but the hydrocarbon chain also becomes more soluble, causing less distortion in the water structure. The balance of these two effects is manifested through a minimum in the cmc. Because there are no explicit solvent molecules in our simulations, this cmc minimum is not properly reproduced. The inclusion of explicit solvent molecules would significantly increase the necessary CPU time. In that case, we would not be able to capture the micellization process and obtain estimates of the cmc with our present computer capabilities. The simulation results for the cmc also appear to be somewhat dependent on system size. The cmc at L ) 79.4 Å is significantly higher than at larger box sizes, ranging from about 13 to 38% higher than the values at L ) 119.1 Å. The system size effect seems to be weaker at larger box size. The cmc values at L ) 119.1 Å are about 6 to 27% higher than the values at L ) 158.8 Å. This suggests that the cmc converges at large system sizes. The cmc appears to be depend roughly linearly on 1/L, although we do not have enough data to extrapolate the cmc confidently to the infinite size limit. To investigate the effect of system size further, we have looked at the dependence of the volume fraction on the chemical potential (data not shown). At a high temperature T ) 500 K, system size has very little effect. At a lower temperature (though still relatively high) of T ) 400 K, the curve for L ) 79.4 Å deviates significantly from the larger box sizes. However, the curves for L ) 119.1, 158.8, and 238.2 Å show good agreement. Therefore, the discrepancy in the cmc for L ) 79.1 Å can probably be attributed to system size effects, and further simulations to investigate the cluster size distribution and counterion binding were performed for L ) 158.8 Å. First-order macroscopic phase transitions also exhibit sharp breaks in the thermodynamic properties such as the osmotic pressure versus concentration curve that we are using to determine the cmc. Therefore, it is important to differentiate between micellization and macroscopic phase separation from simulations at a fixed system size. We have confirmed that our model SDS undergoes micellization by investigating the effect of system size on the apparent phase diagrams, as described in a previous study by our group of nonionic surfactants.56 The apparent volume fraction in the surfactant-rich “phase” shows a major decrease with increasing simulation box size, as expected for systems that undergo micellization.
5. Cluster Size and Counterion Binding We have performed an analysis of the cluster size distribution of the model SDS as a function of temperature and overall surfactant volume fraction. The volume fraction of clusters as a function of the cluster aggregation number M at T ) 350 K for varying overall surfactant volume fractions in L ) 158.8 Å is shown in Figure 6. The cluster size distribution shows two peaks, one at a low aggregation number corresponding to the monomers and another at a higher aggregation number corre-
Langmuir, Vol. 22, No. 9, 2006 4081
Figure 6. Volume fraction of clusters φM vs the cluster aggregation number M at T ) 350 K for various overall surfactant volume fractions in a simulation box of L ) 158.8 Å. Black corresponds to 48 mM, red to 26 mM, blue to 10 mM, and pink to 5 mM.
sponding to the micelles. At this temperature, we see that the aggregation number is somewhat sensitive to the overall surfactant concentration, increasing from M ) 26 at 26 mM to M ) 31 at 48 mM at T ) 350 K. The aggregation number shows a similar increase from M ) 23 to 26 as the concentration goes from around 48 mM to about 76 mM at a higher temperature of T ) 400 K (data not shown). The concentrations are much higher than the cmc because of the relatively small box size. Because the cmc occurs at such low concentrations, very large boxes are necessary for micelles to exist at concentrations around the cmc. From these distributions, we see that even when the surfactant concentration is almost doubled the average aggregation number increases slightly only, suggesting that the addition of surfactants leads to the formation of additional micelles of similar sizes rather than the growth of existing clusters. The aggregation number is also dependent on temperature, increasing from M ) 23 to 31 at roughly the same concentration (48 mM) as the temperature decreases from T ) 400 to 350. At T ) 300 K, micelles of around M ) 40 are formed as shown in Figure 3. These aggregation numbers of M ≈ 20-40 are smaller than the aggregation number of SDS determined experimentally, which is in the range of M ≈ 60-80 at room temperature.8,85 The discrepancy is probably due to the simplification of the hydrocarbon structure. The approximation of two methyl groups as a single lattice site probably could not accurately reproduce the structure of the micellar core. The excluded volume is overestimated when using a lattice model, resulting in the formation of smaller micelles. We have also investigated the dependence of the concentration of free surfactants and of surfactants that exist as micelles on the overall surfactant concentration. The results at T ) 400 K in L ) 158.8 Å are shown in Figure 7. The concentration of free surfactants is defined as the fraction of surfactants that exist as aggregates up to the cluster size at the minimum between the two peaks of the cluster size distribution (Figure 6). The rest of the surfactants are considered to exist in micelles. From Figure 7, we see that the concentration of free surfactants increases with the overall concentration until a certain threshold. It then levels off and decreases slightly at high overall surfactant concentrations. The concentration at which the free surfactant concentration levels off is sometimes considered to be the cmc. This break occurs at around 10 mM, in good agreement with our estimate of the cmc at 400 K. The decrease in free surfactant concentration at high (85) Chen, J. M.; Su, T. M.; Mou, C. Y. J. Phys. Chem. 1986, 90, 2418.
4082 Langmuir, Vol. 22, No. 9, 2006
Figure 7. Concentration of free surfactants and surfactants in micelles as a function of the overall surfactant concentration at T ) 400 K in a simulation box of L ) 158.8 Å. The circles correspond to the free surfactants, and the triangles correspond to surfactants that exist in micelles.
overall surfactant concentration can be attributed to an excluded volume effect. As the concentration of micelles increases, there is less available volume for the free surfactants. This decrease has also been observed in previous experimental86,87 and simulation47,55 studies. This effect has also been incorporated in theories of micelle formation.23 The concentration of surfactants in micelles increases linearly after the cmc is crossed. Another property of interest is the degree of counterion binding to the charged micelles. In this study, a counterion is considered to be bound to a micelle if the distance dij between the counterion and a headgroup of a surfactant that belongs to a cluster of at least 20 surfactants satisfies the constraint σht e dij < (σht + 1). By this convention, a system with no micelles will not have any bound counterions, and counterions associated with free surfactant chains (not part of a micelle) are not included. The definition of a bound counterion is somewhat arbitrary, and the degree of counterion binding will be sensitive to the choice of the cutoff distance. However, it will not affect the qualitative behavior of the counterion binding. The fraction of counterions considered to be bound to micelles was determined for L ) 158.8 Å at different temperatures and overall surfactant volume fraction and is shown in Figure 8. Experimental values for the counterion binding obtained from a potentiometric method88 at T ) 298 K are also shown in Figure 8. The potentiometric method allows the direct measurement of the concentration of free counterions in solution. The potential measured in the surfactant solution is compared with the potential of an electrolyte solution at the same concentration. At low surfactant concentration (below the cmc), no micelles exist, and the two potentials are equivalent. Above the cmc, micelles begin to form, and counterion condensation occurs. There is a difference in the potentials, which is directly related to the degree of counterion binding. This method also allows an estimation of the cmc. Because counterion binding occurs only when micelles are present in the solution, the concentration at which the measured counterion concentration begins to be less than the overall concentration will be approximately the cmc. The cmc was estimated experimentally to be approximately 8 mM by this method,88 in good agreement with cmc values obtained by other methods.10,16,84 Figure 8 shows that above the cmc the degree (86) Johnson, I.; Olofsson, G.; Jo¨nsson, B. J. Chem. Soc., Faraday Trans. 1 1987, 83, 3331. (87) Kahlweit, M.; Teubner, M. AdV. Colloid Interface Sci. 1980, 13, 1. (88) Hsiao, C. C.; Wang, T.-Y.; Tsao, H.-K. J. Chem. Phys. 2005, 122, 144702.
Cheong and Panagiotopoulos
Figure 8. Degree of counterion binding as a function of overall surfactant concentration and temperature. Open circles are experimental values obtained from potentiometric methods by Hsiao et al.88 Filled symbols are the results from this study. Circles, squares, and triangles are for T ) 300, 350, and 400 K, respectively, in L ) 158.8 Å. Open triangles are also from this work for T ) 400 K in L ) 238.2 Å. Lines are guides for the eye.
of counterion binding increases steadily before leveling off at a constant value of about 0.785 at a concentration between 50 and 100 mM. The simulation results are in good qualitative agreement with the experimental results. However, at high concentrations, the degree of counterion binding is lower than the experimental values. The numerical discrepancies can be attributed to a few possible reasons. The degree of binding increases with an increase in the micelle size. A micelle consisting of more surfactant molecules means that the micelle charge is larger. A higher micelle charge is likely to attract more counterions, increasing the degree of counterion binding. As discussed previously, the lattice model used in this study underpredicts the micelle size, so it is expected that the degree of counterion binding is also lower. In addition, the definition of a bound counterion is somewhat arbitrary in our simulations. If the cutoff distance that defines a bound counterion is increased, then the degree of counterion binding will also increase. At low concentrations, counterion condensation occurs at lower concentrations than the experimental values. This is a result of the low predicted critical micelle concentration by the simulations. Counterion condensation occurs only when micelles are present. Because micelles form at lower concentrations in the simulations, the predicted onset of counterion condensation also occurs at lower concentrations. Despite the numerical discrepancies, the qualitative behavior is not expected to change. Counterion binding is found to depend on temperature. At low temperatures, the electrostatic interactions are strong, increasing the attractions between counterions and micelles. Thus, we find a higher degree of counterion binding at T ) 300 K. At higher temperatures, the electrostatic interactions are weaker, resulting in a lower degree of counterion binding. This dependence on temperature has also been observed in experimental studies.89,90 We have obtained the counterion binding at T ) 400 K for two system sizes, L ) 158.8 and 238.2 Å. The results were found to be within the simulation uncertainties, so system size effects are not significant for counterion condensation. (89) Oremusova, J.; Mislovicova, V. Collect. Czech. Chem. Commun. 1997, 62, 1853. (90) Moroi, Y.; Take’uchi, M.; Yoshida, N.; Yamauchi, A. J. Colloid Interface Sci. 1998, 197, 221.
Micellization in Model Ionic Surfactants
Langmuir, Vol. 22, No. 9, 2006 4083
6. Conclusions
poulos,91 TraPPE,92 NERD,93 and OPLS,94 have been developed and optimized to the various thermodynamic properties of alkanes. It is possible to combine this model for dodecane in continuum space, with a finely discretized lattice model of the headgroups and counterions. The lattice will allow the long-range Coulombic interactions to be precomputed, whereas the realistic hydrocarbon model will allow the tail structure to be better reproduced. It would also be of great interest to perform a systematic study on the effect of the various parameters on the micellization properties. The size and valence of the headgroup, length of the tail, dielectric constant, and strength of interactions can all be systematically varied to understand their effect and importance in micellization. Other interesting variations include allowing for multivalent counterions and varying the position of the charged group along the linear chain. These factors are expected to play a large role in the micellization properties of a surfactant. Also of great interest would be the micellization behavior of surfactant solutions in the presence of salt. Investigating the phase behavior of mixtures of surfactants (i.e., ionic and nonionic surfactant mixtures) would also have great practical applications. An eventual goal of these simulations would be the inclusion of explicit solvent molecules and the realistic simulations of surfactant and micellar solutions.
We have made progress toward developing a model for ionic surfactants that allows us to predict their micellization properties quantitatively. The strength of the model lies in its simplicity and generality, and it can be used to study different surfactants. We obtained parameters for sodium dodecyl sulfate by matching the solubility of hydrocarbons in water at 300 K and the liquid density of dodecane. It was found that by allowing a tail bead to represent two carbon groups and choosing an appropriate energy parameter for the tail-tail interactions and bead diameter properties such as the solubility of hydrocarbons in water, the liquid density, the end-to-end distance, and the radius of gyration, the parameters are reasonably matched. Using this set of parameters, we performed grand canonical Monte Carlo simulations to obtain the micellization properties of SDS. The properties of interest are the cmc, the average aggregation number, and the degree of counterion binding. Predictions of the various micellization properties exhibit the correct qualitative trends and give quantitative values in the right order of magnitude when compared to experimental values. However, because of its simplicity, the model has several shortcomings. The cmc, micelle size, and degree of counterion binding are all significantly underpredicted, and certain features of the dependence of the cmc on temperature are overlooked by the model. The model can be improved to give better agreement with experimental results. The first improvement to the model could be the inclusion of short-range head-head and head-tail repulsions. Because head-head repulsions raise the cmc, this could potentially alleviate the problem of underpredicting the cmc. It would also be of great practical interest to develop a more realistic model for ionic surfactants that can quantitatively predict the micellization properties of a range of surfactants. For the model to be more realistic, it is important to represent the hydrocarbon chain as well as possible. A possible strategy would be to use a more detailed description of the hydrocarbon chain. Numerous potentials, such as that of Errington and Panagioto-
Acknowledgment. Funding for this work was provided by the Department of Energy, Office of Basic Energy Sciences (DEFG201ER15121). Additional support was provided by ACSPRF (grant no. 38165-AC9). Note Added after ASAP Publication. This article was released ASAP on March 9, 2006. A y-axis label was added to Figure 8, and the correct version was posted on March 20, 2006. LA053511D (91) Errington, J. R.; Panagiotopoulos, A. Z. J. Chem. Phys. 1999, 103, 6314. (92) Martin, M. G.; Siepmann, J. I. J. Chem. Phys. 1998, 102, 2569. (93) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1998, 108, 9905. (94) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225.
