746
Langmuir 1996, 12, 746-753
Simulations of Micelle Self-Assembly in Surfactant Solutions Bruce J. Palmer* and Jun Liu Energy and Environmental Sciences Division, Pacific Northwest National Laboratory, Richland, Washington 99352 Received July 7, 1995. In Final Form: October 24, 1995X A simple model of surfactants in aqueous solution is described and molecular dynamics simulations are presented. This model qualitatively reproduces the energy ordering of interactions for surfactants in solution, instead of using the symmetric solvent-solvent and hydrocarbon-hydrocarbon interactions used previously in similar models. The effects of surfactant chain length and concentration are investigated. The results show that the tendency toward micelle formation in the model increases as the chain length is increased and that the type of micelle formed varies with surfactant concentration. Transitions from loose aggregates to spherical micelles to cylindrical micelles are seen as surfactant chain length is increased. A transition from cylindrical micelles to spherical micelles is also seen as surfactant concentration is lowered.
1. Introduction The phenomenon of self-assembly in surfactant solutions plays an important role in biological and industrial processes such as food digestion, chemical separation, and surface catalysis.1,2 Recently, surfactant molecules have been used to mediate and direct the synthesis of inorganic materials with well-defined microstructures and morphologies.3,4 Novel zeolite-like mesoporous molecular sieves have been synthesized based on hexagonal and cubic surfactant aggregates. These materials are characterized by ordered porosity, narrow pore size distribution, and adjustable pore size (from 1.5 to 20 nm).5-7 To understand surfactant self-assembly in such complex solutions, it is desirable to establish theoretical models that realistically account for the major inter- and intramolecular interactions in the surfactant solution. Although surfactant solutions are playing an increasingly important role in many chemical processes, theoretical understanding and modeling of their behavior are still fairly primitive. Over the last decade or so, many researchers have investigated lattice models of emulsions, and there has been considerable progress in mapping out the qualitative features of the phase diagrams for these systems.8-12 More recently, workers have begun looking at continuum models of emulsions and surfactant solutions. Telo da Gamma and Gubbins have used liquid state density functional theory to calculate the surface tension, absorption, and orientation profile for a liquid-liquid X Abstract published in Advance ACS Abstracts, January 15, 1996.
(1) Israelachvili, J. Intermolecular and Surface Forces, 2nd ed.; Academic Press: San Diego, CA, 1991; Chapter 17. (2) Evans, D. F.; Wennerstrom, H. The Colloidal Domain Where Physics, Chemistry, Biology, and Technology Meet; VCH Publishers, Inc.: New York, 1994; Chapter 1. (3) Barnickel, P.; Wokaun, A. Mol. Phys. 1990, 69, 1. (4) Lisiecki, I.; Pileni, M. P. J. Am. Chem. Soc. 1993, 115, 3887. (5) Beck, J. S.; Vartuli, J. C.; Roth, W. J.; Leonowicz, M. E.; Kresge, C. T.; Schmitt, K. D.; Chu, C. T.-W.; Olson, D. H.; Sheppard, E. W.; McCullen, S. B.; Higgins, J. B.; Schlenker, J. L. J. Am. Chem. Soc. 1992, 114, 10834. (6) Kresge, C. T.; Leonowicz, M. E.; Roth, W. J.; Vartuli, J. C.; Beck, J. S. Nature 1992, 359, 710. (7) Monnier, A.; Schuth, F.; Huo, Q.; Kumar, D.; Margolese, D.; Maxwell, R. S.; Stucky, G. D.; Krishnamurty, M.; Petroff, P.; Firouzi, A.; Janicke, M.; Chmelka, B. F. Science 1993, 261, 1299. (8) Wheeler, J. C.; Widom, B. J. Am. Chem. Soc. 1968, 90, 3064. (9) Schick, M.; Shih, W.-H. Phys. Rev. B 1986, 34, 1797. (10) Halley, J. W.; Kolan, A. J. J. Chem. Phys. 1988, 88, 3313. (11) Widom, B. J. Chem. Phys. 1984, 81, 1030. (12) Widom, B. J. Chem. Phys. 1986, 84, 6943.
0743-7463/96/2412-0746$12.00/0
interface in a three-component system composed of two immiscible liquids and a bipolar amphiphile.13 The continuum studies have been extended to slightly more detailed models of surfactant solutions and emulsions using computer simulations to examine micelle formation and the behavior of surfactants at the oil-water interface.14-19 Rector et al. have been able to map out the thermodynamics of cluster formation using simulations of this type on a simple two-site surfactant.20 The surfactants in these models are composed of chains of Lennard-Jones sites connected by harmonic springs, and the solvent is modeled as a simple Lennard-Jones fluid. These models treat the hydrocarbon-hydrocarbon and water-water interactions symmetrically, in the sense that the Lennard-Jones potentials for these two interactions are the same. The hydrocarbon-water interaction is modeled using a short-range repulsive potential. The rationale for making the hydrocarbon-hydrocarbon interaction as large as the water-water interaction has been that the hydrocarbon sites in these models are meant to represent more than one CH2 group. Although these models still represent a crude picture of surfactants in aqueous solution, they remain attractive candidates for study because of the ease with which large systems can be simulated over long time periods with relatively modest computer resources. Simulations have also been performed using detailed molecular mechanics potentials designed to realistically model the behavior of actual surfactant solutions. Watanabe et al. have done simulations on solutions of sodium octanoate in water21-23 and Egberts and Berendsen have reported results for a sodium decanoate/decanol/water (13) Telo da Gama, M. M.; Gubbins, K. E. Mol. Phys. 1986, 59, 227. (14) Smit, B. PhD Thesis, Rijksuniversiteit Utrecht, The Netherlands, 1990. (15) Smit, B.; Hilbers, P. A. J.; Esselink, K.; Rupert, L. A. M.; van Os, N. M.; Schlijper, A. G. Nature 1990, 348, 624. (16) Smit, B.; Hilbers, P. A. J.; Esselink, K.; Rupert, L. A. M.; van Os, N. M.; Schlijper, A. G. J. Phys. Chem. 1991, 95, 6361. (17) Karaborni, S.; van Os, N. M.; Esselink, K.; Hilbers, P. A. J. Langmuir 1993, 9, 1175. (18) Smit, B.; Esselink, K.; Hilbers, P. A. J.; van Os, N. M.; Rupert, L. A. M.; Szleifer, I. Langmuir 1993, 9, 9. (19) Karaborni, S.; Esselink, K.; Hilbers, P. A. J.; Smit, B.; Karthauser, J.; van Os, N. M.; Zana, R. Science 1994, 266, 254. (20) Rector, D. R.; van Swol, F.; Henderson, J. R. Mol. Phys. 1994, 82, 1009. (21) Watanabe, K.; Ferrario, M.; Klein, M. L. J. Phys. Chem. 1988, 92, 819. (22) Watanabe, K.; Klein, M. L. J. Phys. Chem. 1989, 93, 6897. (23) Watanabe, K.; Klein, M. L. J. Phys. Chem. 1991, 95, 4158.
© 1996 American Chemical Society
Micelle Self-Assembly
Langmuir, Vol. 12, No. 3, 1996 747
system.24 Alper et al. have also done simulations on phospholipid monolayers and membranes.25,26 While these simulations provide a great deal of detailed information about real systems, they require large investments in computer time in order to run even a single simulation. This paper will take the approach of using quasi-realistic models of surfactants built up from a combination of Lennard-Jones sites and harmonic bond-stretching and angle-bending potentials. These models are easy to study numerically and allow multiple simulations to be performed under a variety of conditions. The most significant difference between the models investigated here and those used by previous authors is that the symmetric treatment of the hydrocarbon-hydrocarbon and water-water interactions is dropped and a more realistic ordering of the relative energetics of these interactions is used instead. Preliminary investigations are performed on the surfactant chain length and effect of concentration on micelle formation. 2. The Model The model surfactant systems investigated in this paper consist of simple surfactants dissolved in a Lennard-Jones solvent. The surfactant molecules are composed of a head site attached to a chain of tail sites, with the head sites and tail sites bonded to each other via simple harmonic springs. The different sites along the surfactant chain are also subject to a harmonic bending potential that favors 180° angles in the tail. Using a 180° equilibrium angle instead of 120° increases the effective length of the surfactant by preventing any sharp bends from occurring in the surfactant tail, but it also prevents the excessive folding that occurs if no bending potential is used. Simulations by Watanabe et al.21,23 and Egberts and Berendsen24 indicate that 70-80% of the dihedral angles in the hydrocarbon tail of surfactants in micelles are in the trans configuration, which supports the use of a straight tail. The sites in the surfactant interact with sites in other surfactant molecules and in the solvent via Lennard-Jones interactions. The potential function for the model surfactant solution is a standard molecular mechanics potential of the form N
U)
φij(rij) + ∑ K(θijk - π)2 + ∑ B(rij - R0)2 ∑ i Rijc (2.2)
The term φ0(Rcij) is a constant chosen so that the pair interaction vanishes at the cutoff potential distance Rcij. Two types of cutoff are used in these calculations. (24) Egberts, E.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 3718. (25) Alper, H. E.; Bassolino-Kima, D.; Stouch, T. R. J. Chem. Phys. 1993, 99, 5547. (26) Alper, H. E.; Stouch, T. R. J. Phys. Chem. 1995, 99, 5724.
Table 1. Lennard-Jones Interaction Parameters for the Surfactant-Solvent System interaction
Rijc /σij
σij
ij
solvent-solvent solvent-head solvent-tail head-head head-tail tail-tail
2.5 2.5 2.5 21/6 2.5 2.5
1.0 1.5 1.0 2.0 1.5 1.0
1.0 1.0 0.125 1.0 0.125 0.125
Table 2. Lennard-Jones Interaction Parameters for the Surfactant-Solvent System with Strong Tail-Tail Interactions interaction
Rijc /σij
σij
ij
solvent-solvent solvent-head solvent-tail head-head head-tail tail-tail
2.5 2.5 21/6 21/6 21/6 2.5
1.0 1.5 1.0 2.0 1.5 1.0
1.0 1.0 1.0 1.0 1.0 1.0
The first is Rcij ) 2.5σij and results in a conventional cutand-shifted Lennard-Jones potential with a potential minimum of approximately ij and a hard-sphere radius of σij. The second cutoff is Rcij ) 21/6σij and results in a short-range, purely repulsive interaction. The parameters used in the Lennard-Jones interactions were chosen to qualitatively mimic the behavior of surfactants in aqueous solution. All parameters are defined in terms of the solvent-solvent interaction parameters, which are set equal to σSS ) 1.0 and SS ) 1.0. The solvent-solvent interactions are attractive with the c /σSS ) 2.5. All other Lennard-Jones cutoff set at RSS interactions are defined relative to these numbers and are summarized in Table 1. The hard-sphere radius for the head group is chosen to be twice as large as the sites in either the tail or the solvent. The interaction energy between the head and the solvent is comparable to the solvent-solvent interaction energy, corresponding to a polar or charged head group in aqueous solution. The tail sites interact much more weakly with other sites, corresponding to the dispersion interactions that would be expected for hydrocarbons interacting either with themselves or with other polar species. The choice of a weakly interacting tail is the most significant difference between the models investigated here and those used by Smit et al. to investigate the behavior of surfactant systems.14-19 Smit et al. chose parameters such that the tail-tail interaction and solvent-solvent interactions were of comparable strength and then made the tail-solvent interaction repulsive. These models form micelles quite readily under a variety of conditions and have been used to model both micelles and oil-water emulsions. However, there may be some significant qualitative differences between these strong tail interaction models and the models investigated here. For comparison with the strong tail interaction models,14-19 some simulations were performed using a set of interaction parameters with a strong tail-tail interaction and a repulsive tail-solvent interaction. This system is not exactly the same as those investigated by other workers but is qualitatively very similar. The parameters are summarized in Table 2. The remaining parameters for the bond angle bends, bond stretches, and masses are listed in Table 3. The bond stretching parameter B is strong enough so that the root mean square fluctuation in the bond length is about 5% of the average value. The trajectories are integrated using the velocity Verlet algorithm recast as a three-point Gear predictor-corrector.27,28 The time step was set at 0.005τ0, where τ0 ) (mSσ2SS/SS)1/2. The nonbonded inter-
748
Langmuir, Vol. 12, No. 3, 1996
Palmer and Liu
Figure 1. Configuration from simulation of 100 S8 surfactant molecules in 4000 solvent molecules: (a, top) with solvent; (b, bottom) without solvent. The head sites are red, the tail sites are blue, and the solvent sites are yellow.
actions where calculated using a linked-cell algorithm to construct a neighbor list.29 Most of the simulations were done under constant pressure-constant temperature
conditions, using the Andersen-Nose´ extended system formalism.30,31 To minimize finite size effects, periodic boundary conditions were employed in all simulations.
Micelle Self-Assembly
Langmuir, Vol. 12, No. 3, 1996 749
Figure 2. Configuration from simulation of 100 S4 surfactant molecules in 4000 solvent molecules. The head sites are red and the tail sites are blue. Four copies of the unit cell are displayed. Table 3. Parameters for Intramolecular Potentials and Masses parameter
value
parameter
value
K B R0
10.0 200.0 1.0
mS mH mT
1.0 2.0 1.0
The temperature for all runs was set at T/kBSS ) 0.8 and 3 the pressure was set at P/(SS/σSS ) ) 0.1. This corresponds to a temperature below the critical point and a pressure above the vapor pressure.14 The values of the volume mass MV and the temperature scaling mass MT in the extended system Hamiltonian were set at 0.01 and 1000.0, respectively. Using these parameters, the ratio of the root mean square fluctuations in the Hamiltonian to the average value of the kinetic energy was less than 1 part in 10 000. For this algorithm, such conservation is a good indication that the pressure calculation is being performed correctly. (27) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (28) Palmer, B. J.; Garrett, B. C. J. Chem. Phys. 1993, 98, 4047. (29) Auerbach, D. J.; Paul, W.; Bakker, A. F.; Lutz, C.; Rudge, W. E.; Abraham, F. F. J. Phys. Chem. 1987, 91, 4881. (30) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (31) Nose´, S. J. Chem. Phys. 1984, 81, 511.
3. Results Simulations were initially performed on solutions containing surfactants with four, six, and eight tail sites. These surfactants are referred to in the following as S4, S6, and S8, respectively. A starting configuration containing 100 identical surfactant molecules and 4000 solvent molecules was prepared with the surfactant molecules randomly distributed throughout the system. This was accomplished by starting the system at very low density with the surfactant and solvent molecules randomly distributed on a lattice and then letting the system condense to a liquid solution using constant pressureconstant temperature molecular dynamics. The initial configuration was taken shortly after the system had condensed, when the distribution of surfactant molecules was still fairly random. This configuration was then allowed to evolve under constant pressure-constant temperature conditions in blocks of 50 000 steps. A configuration from the simulation of the S8 surfactant is shown in Figure 1, where the same configuration is displayed both with and without the solvent. The solvent generally tends to obscure the details of the surfactant structure, so it is not included in any of the remaining figures. Figure 1 shows the formation of a cylindrical micelle with the view straight down the axis of the micelle.
750
Langmuir, Vol. 12, No. 3, 1996
Palmer and Liu
Figure 3. Configuration from simulation of 100 S6 surfactant molecules in 4000 solvent molecules. The head sites are red and the tail sites are blue. Four copies of the unit cell are displayed.
From the figure, it is not immediately clear that a single micelle has formed because the center of the cylinder happens to lie on one of the boundaries of the periodic cell. To make surfactant structures that appear at the faces and corners of the simulation cell easier to visualize, four replicas of the fundamental simulation cell are shown in the remaining figures. While this makes the identification of surfactant structures easier, it also introduces an artificial lattice structure to the configuration which should be ignored. Configurations from the S4, S6, and S8 surfactant simulations are shown in Figures 2-4. All simulations contain 100 surfactant molecules and 4000 solvent molecules. The configurations represent systems that were run for at least 1.5 × 105 time steps. The configuration for the S4 surfactant shows some signs of weak aggregation, but the structures appear fairly loosely knit and there are no clear indications of a recognizable micelle forming. The S6 surfactant configuration, however, appears to contain at least one large globular aggregate, along with a smaller satellite aggregate. Because the graphical visualization projects everything into the plane and gives a poor sense of depth, all configurations were checked by viewing along the three coordinate axes to make sure that the apparent structures were not artifacts
of the visualization process. The two clusters visible in Figure 3 showed up in all three views. Finally, Figure 4 shows two views of the cylindrical micelle formed in the simulation of the S8 surfactant system. Figure 4a is a side view of the micelle, showing that it is continuous across the simulation cell boundaries. Figure 4b is equivalent to the view down the micelle axis shown in Figure 1, except that it includes four replicas of the simulation cell. A ring of head groups is clearly visible surrounding a core composed to tail sites. This micelle looks remarkably like the cylindrical sodium octanoate micelle studied by Watanabe and Klein.23 From Figure 1a it can also be seen that the solvent is almost completely excluded from the micelle core. A question that immediately arose in regard to the formation of the cylindrical micelle in the simulation of the S8 surfactant is whether or not this system would produce spherical micelles if the concentration of surfactant was decreased. To answer this, two additional simulations were performed in which the surfactant concentration was lowered by a factor of 2. In the first simulation, the concentration was lowered by doubling the number of solvent molecules to 8000. In the second, the number of solvent molecules remained fixed at 4000 and the number of surfactant molecules was halved to 50.
Micelle Self-Assembly
Langmuir, Vol. 12, No. 3, 1996 751
Figure 4. Configuration from simulation of 100 S8 surfactant molecules in 4000 solvent molecules: (a, top) side view of cylindrical micelle; (b, bottom) view down the axis of cylindrical micelle. The head sites are red and the tail sites are blue. Four copies of the unit cell are displayed.
Both simulations resulted in the formation of spherical micelles. The simulation containing 8000 solvent mol-
ecules produced two spherical micelles of approximately equal size, and the simulation containing 50 surfactant
752
Langmuir, Vol. 12, No. 3, 1996
Palmer and Liu
Figure 5. Configuration from simulation of 100 S4′ surfactant molecules in 4000 solvent molecules. The head sites are red and the tail sites are blue. Four copies of the unit cell are displayed.
molecules resulted in the formation of a single large micelle containing most of the surfactant molecules in the system. Although these simulations are suggestive of the presence of a phase boundary between cylindrical and spherical phases, they do not constitute conclusive proof of a phase transition. Determining the presence of such a phase boundary would require considerably more simulations. In addition to the transition from cylindrical to spherical micelles, it would be interesting to establish whether there is a plateau in the micelle size as the surfactant concentration is lowered further. Such a plateau is seen experimentally for many systems as the concentration is lowered toward the critical micelle concentration. However, determining whether such a plateau exists would be significantly more challenging than the simulations presented here, both because of the larger number of solvent sites required and because the relaxation times would increase as the surfactant concentration became more dilute. To compare the simulations performed here with the work of other groups, a simulation of the surfactant system with the strongly interacting tail sites was also performed, using a surfactant with four tail sites (this surfactant will be referred to as S4′). A total of 100 S4′ surfactants and 4000 solvent molecules were used in the simulation. After 105 time steps, the system appeared as shown in Figure
5. There are four distinct micelles present in the system, even though the corresponding simulation using the S4 surfactant showed only loosely structured aggregates. This is in agreement with the simulations of Smit et al., who obtained well-defined structures even at fairly high temperatures.17-19 Besides identifying micelle formation via direct visual examination using graphics software, statistics on cluster formation were also collected. Two criteria were investigated for determining whether a pair of surfactant molecules belonged to the same cluster. The first was an energy criterion: two surfactants were considered to be in the same cluster if the total tail-tail interaction energy was less than some energy Ecl. The second was a distance criterion: two surfactant molecules were considered to be in the same cluster if any pair of sites on the two surfactant tails were within some distance Rcl of each other. Because the simulation shown in Figure 4 produced a single micelle that appeared to comprise all the surfactants in the system, it was used as the basis for determining the values of Ecl and Rcl. The values of Ecl and Rcl were chosen so that any further restriction on the pair criteria resulted in a cluster with less than 100 surfactant molecules for the cylindrical micelle system. Using this approach, the values of Ecl ) -0.01SS and Rcl ) 2.0σSS were finally chosen. Cluster distributions for the simulations of 100 sur-
Micelle Self-Assembly
Langmuir, Vol. 12, No. 3, 1996 753
that fall outside the pair critiria for being in the same cluster. Similar behavior was also seen for the cluster distribution calculated from the simulation on the S8 surfactant. This distribution had a sharp peak centered at 100 molecules, corresponding to the single cylindrical micelle seen in Figure 4. 4. Conclusions
Figure 6. Cluster distribution statistics from simulation of 100 S4 surfactant molecules in 4000 solvent molecules: (solid line) energy criteria; (dotted line) distance criteria.
Figure 7. Cluster distribution statistics from simulation of 100 S6 surfactant molecules in 4000 solvent molecules: (solid line) energy criteria; (dotted line) distance criteria.
factants in 4000 solvent molecules are shown in Figures 6 and 7 for the S4 and S6 simulations. The x axis depicts the size of the cluster while the y axis lists the average number of molecules found in a cluster of size x. The distributions were calculated over 50 000 time steps. For both simulations, the energy and distance criteria for two surfactant molecules to be in the same cluster give almost identical distributions. The S4 distribution in Figure 6 shows three distinct peaks centered at 18, 27, and 35 molecules, as well as a significant number of free monomers. However, the presence of a broad background in the distribution and the width of the peaks suggests that the clusters are partially breaking up and re-forming during the course of the simulation. This behavior would be consistent with the loosely structured aggregates seen in Figure 2. The cluster distribution for the S6 simulation is shown in Figure 7 and shows two sharp peaks at 16 and 84 molecules, as well as a peak at 100 molecules. The distribution indicates that there are almost no free monomers in the system. The two peaks are completely consistent with the two aggregates seen in Figure 3. The peak at 100 molecules suggests that the two aggregates are occasionally bumping into each other or may even be starting to merge to form a single micelle (the configuration in Figure 3 was taken at the beginning of the run during which the distribution was calculated). The low background level of other size clusters indicates that the micelles are temporarily breaking up into separate pieces
The simulations described here are clearly capable of producing micelles using interactions that qualitatively reproduce the energy ordering of interactions expected in aqueous solutions. The tendency toward micelle formation increases as the chain length of the surfactant is increased, which reproduces the experimental trend.1 The simulations also suggest that the type of micelle formed changes as a function of chain length and surfactant concentration. The model is sophisticated enough to simulate selfassembly in more complicated systems. The comparison of the simulations on the S4 and S4′ surfactant systems indicates that the micelles formed in the simulations described in this paper are much more loosely organized than the micelles seen in simulations by previous researchers using a strong tail-tail interaction. Because of their stronger binding, the surfactants with strong tail-tail interactions may overestimate the tendency to form micelles, and the micelles they do form may be much more influenced by geometric packing considerations than by such thermodynamic conditions as concentration and temperature. The sharp peaks visible in all the cluster distributions indicate that the micelles formed in these simulations are stable on a time scale comparable to 50 000 steps. The long lifetime of the clusters, combined with the small number of clusters present in each of the simulations, indicates that the calculated distributions do not reflect the equilibrium distribution that would be seen for an infinite system. To calculate distributions that accurately reflect the range of cluster sizes will probably require simulations on much larger systems and may need to be combined with heating-quenching cycles in order to increase the interconversion rates between different size micelles. The simulations, particularly on the S8 surfactant system, point to the need for additional diagnostics to determine whether or not micelles have formed, as well as methods for analyzing micelle shape. Aside from cluster distributions, most current methods for analyzing micelles in simulations appear to rely on graphics visualization to first identify a micelle. The molecules in this micelle can then be tagged and their subsequent behavior during the simulation can be tracked. Future efforts will be devoted to developing diagnostics that do not require prior visual identification and that can be related to experimental probes of micelle formation. The models investigated here will also be extended to include solutes, and the effects of solutes on self-assembly will be investigated. Acknowledgment. The authors are grateful to Dave Rector for many useful discussions and to Fadel Erian for his encouragement and support. This work was funded under the Laboratory Directed Research and Development program as part of Pacific Northwest National Laboratory’s Advanced Processing Technology Initiative. Pacific Northwest Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06-76RLO 1830. LA950979F
