J . Phys. Chem. 1984, 88, 6363-6366
6363
Internal Structure of a Model Micelle via Computer Simulation J. M. Haile* Department of Chemical Engineering, Clemson University, Clemson, South Carolina 29631
and J. P. O’Connell Department of Chemical Engineering, University of Florida, Gainesville, Florida 3261 1 (Received: June I I , 1984; In Final Form: August 16, 1984)
The molecular dynamics method has been used in a preliminary computer simulation study of intramicellar structure. The simulation was performed on a simplified model of a spherical micelle composed of 40 13-member chains with head groups fixed on a spherical shell of 16 A radius. Results are reported for the local density of CH2 and CH, groups, for the distribution of gauche bonds, and for the average conformation of the chains.
Introduction Despite intense study, the molecular structure of micellar systems is still controversial. Theoretical models of such structures include radially aligned, fully extended chains,’Y2 bent matchstick^",^ and models involving several gauche bonds per chain.&’ Experimental probes of micellar structure are provided by NMR,8,9 Raman light scattering,1° and small-angle neutron diffra~tion.1~-’~In addition, Monte C a r l 0 ~ ~ and 3 ’ ~ molecular dynamicsI6J7 methods for simulating matter on digital computers are particularly suited for exploring microscopic structure. Apparently, however, other than a Monte Carlo study on a simplified lattice model,18these simulation methods have not been applied to micelles. This paper describes a first attempt to study intramicellar structure using the molecular dynamics method. The micelle model used here, while quite approximate, reveals important features of chain conformations and provides a basis for future, more sophisticated simulation studies. Description of the Model We constructed our first micelle model from 40 skeletal alkane (1) See e.g. C. R. Cantor and P. R. Schimmel, ‘Biophysical Chemistry”, Vol. 111, W. H. Freeman, San Francisco, 1980, p 1346; P. W. Atkins, “Physical Chemistry”, 2nd ed., W. H. Freeman, San Francisco, 1982, p 845; I. N. Levine, “Physical Chemistry”, McGraw-Hill, New York, 1978, p 342. (2) F. M. Menger and P. C. Vesquez, J. Org. Chem., 47,5400 (1982), and references cited therein. (3) P. Fromherz, Chem. Phys. Lett., 77, 460 (15181). (4) K. A. Dill and P. J. Flory, Proc. Natl. Acad. Sei. U.S.A., 77, 3115 (1980); 78, 676 (1981). (5) K. A. Dill, D. E. Koppel, R. S. Cantor, J. A. Dill, D. Bendedouch, and S.-H. Chen, Nature (London), 309, 42 (1984). (6) K. A. Dill in ‘Surfactants in Solution“, K. L. Mittal, et al., Eds., Plenum Press, New York, 1983. (7) D. W. R. Gruen in “Surfactants in Solution”, K. L. Mittal, et al., Eds., Plenum Press, New York, 1983; J . Colloid Interface Sci., 84, 281 (1981). (8) B. J. Cabane, J . Phys. (Orsay, Fr.), 42, 847 (1981). (9) P. G. Nilsson, H. Wennerstrom, and B. Lindman, J . Phys. Chem., 87, 1377 (1983); B. Lindman and H. Wennerstrom in “Topics in Current Chemistry: Micelles”, Vol. 7, Springer-Verlag. New York, 1980. (10) K. Kalyanasundaram and J. K. Thomas, J. Phys. Chem., 80, 1462 (1976). (1 1) D. Bendedouch, S. H. Chen, and W. C. Koehler, J . Phys. Chem., 87, 153 (1983). (12) D. Bendedouch, S. H. Chen, W. C. Koehler, and J. S. Lin, J . Chem. Phys. 76, 5022 (1982). (13) J. B. Hayter and J. Penfold, J. Chem. Soc., Faraday Trans. 1, 77, 1851 (1981); J. B. Hayter and T. Zemb, Chem. Phys. Lett., 93, 91 (1982). (14) K. Binder, Ed., “Monte Carlo Methods in Statistical Physics”, Springer-Verlag, New York, 1979. (IS) J. P. Valleau and S. G. Whittington in “Statistical Mechanics”, Part A, B. J. Berne, Ed., Plenum Press, New York, 1977. (16) J. Kushick and B. J. Berne in “Statistical Mechanics”, Part B, B. J. Berne, Ed., Plenum Press, New York, 1977. (17) F. Vesely, “Computerexperimente an Flassigkeitsmo4eIlenn, Physik Verlag, Vienna, 1978. (18) S. W. Haan and L. R. Pratt, Chem. Phys. Lett., 79,436 (1981); 81, 386, (1981).
0022-3654/84/2088-6363$01.50/0 , I
,
chains, each composed of 13 soft spheres, mimicking a dodecyl surfactant. The soft spheres represent methyl and methylene groups that interact with intermolecular forces through a Lennard-Jones (6,9) pair potential and interact with intramolecular forces through harmonic potentials for bond vibration and bond angle bending. The functional forms and values of potential parameters for these inter- and intramolecular forces were all taken from Weber;I9 in particular, the carbon-carbon bond length in this model is 1.539 A and u, the pair separation at zero potential, is 3.5 A. For the bond rotational potential we used a model due to Ryckaert and Bellemans,20*21 since we believe it to be more realistic than Weber’s rotational potential. Because this was a preliminary study, we did not use explicit head-head and head-solvent interactions; instead, the centers of the 40 head groups were fixed at points within segments of equal area distributed over a spherical surface. The final radius of the sphere was 16 A, which compares with 17.1 A for the sum of the equilibrium, all-trans, center-of-head to center-of-tail distance (15.36 A) plus the tail methyl group radius (1.75 A). One measure of the chain density is the fraction of the total volume Yoccupied by the repulsive portions of the groups. This volume fraction is given by Nm3/6Y, where N is the total number of groups in the model micelle, and in our case its value is 0.65. This value is comparable to volume fractions used in other molecular dynamics studies of chain molecule systems. In the liquid butane work of Ryckaert and Bellemans20~21 this fraction is 0.76, primarily because their u value is larger (3.92 A). In Weber’s butane simulation^,^^ the fraction varies from 0.40 to 0.68. The simulation study of a bilayer, performed by van der Ploeg and Berendsen:2 can be compared only on a two-dimensional basis; their fractional area occupied was 0.44. The heads in our study were allowed to rotate, and soft-sphere interactions were included between head groups and neighboring methylene groups. The solvaphobic effect acting between methylene groups and the solvent was modeled by an r-I2 repulsive shell displaced outward about one CH2-group diameter from the head-group sphere. This displacement allowed chain segments to lie adjacent to or to penetrate through the head-group sphere. Except for the fixed head groups, the model is (coincidentally) similar to that of van der Ploeg and Berendsen.22
Description of the Simulation Newton’s second-order differential equations of motion were solved for each of the 520 soft spheres by using a fifth-order, (19) T. A. Weber, J . Chem. Phys., 69, 2347 (1978). (20) J. P. Ryckaert, D. Sc. Thesis, Libre Universit6 de Bruxelles, 1976. English translation by J. M. Haile and P. W. Haile, University Microfilms Intl., Ann Arbor, MI, 1980. (21) J. P. Ryckaert and A. Bellemans, Chem. Phys. Lett., 30, 123 (1975). (22) P. van der Ploeg and H. J. C. Berendsen, J . Chem. Phys., 76, 3271 (1982).
0 1984 American Chemical Society
site on the surface of a sphere. (2) A momentum scaling procedurez4was used to maintain the total kinetic energy at a constant
/r“\I\\, I3
I/ 3-
I’:
I
va1ue.”’21 After executing’’2’ time step’ the number conformation had decreased from ‘O0% to Of bonds in the 56%* Throughout*a trans bond is defined a‘ One whose angle 4 satisfies cos I$ < -0.50; all others are in a gauche conformation. Over the next 8000 time steps the density of the aggregate was increased by shrinking the micelle volume. At intervals of 10 time steps the positions of the CH2/CH3groups, the head groups, and the repulsive shell were contracted by 0.02 A along a micelle radius. This slow rate of collapse allowed individual chains to arrange themselves within the available volume without assuming highly improbable conformations. When the radius reached 16 A, the collapse was ended and an additional 3000 time steps were run to allow the system to relax from the imposed transient. At this point 13 000 time steps had been run and the fraction of trans bonds was stable at a mean value of 47.5% with fluctuations of f2%. At time step 13 000 the value of the rotational force constant was reset to Ryckaert’s value,2°*21 causing the instantaneous trans fraction to increase from 47% to 7 1% over the next 13 000 time steps. At the final density of the model micelle, the standard molecular dynamics algorithm executed 45 time steps per CPU minute on an IBM 370/3081 computer. At time step 23 000 the execution speed was improved by changing the algorithm to a multiple time-step procedurez5which executed 100 time steps per CPU minute. After a total of 26 000 time steps, corresponding to 52 ps of real time, we judged the system to be equilibrated and an equilibrium run of 62 000 additional time steps (124 ps) was begun. Over the equilibrium segment the average fraction of trans bonds was 71% with fluctuations generally less than *2%. The major exception occurred during the period 84-104 ps into the equilibrium run, in which the instantaneous trans fraction experienced a large fluctuation from 73% to 65% and then back to about 70%, where it remained as in the earlier part of the equilibrium calculation. Ommended
Results The equilibrium portion of the simulation was analyzed to answer three basic questions about the structure of the model micelle. (1) What is the average radial distribution of methylene and tail groups? (2) What is the average distribution of gauche bonds? (3) What is the variation in a chain’s conformation relative to the micelle radius through the chain’s head group? To answer the first question, the singlet distribution function pi(r) for each CH2/CH3 group was determined. Here i = 2, 3, ..., 13, where group 2 is immediately adjacent to the fixed head (i = 1) and group 13 is the methyl tail. The function pl(r) is proportional to the probability of finding group i at a distance r from the center of the micelle. We have normalized pi(r) so
where dr = 4=+ dr is a spherical volume element and N (= 40) is the total number of chains in the system. The for each (23) C. W. Gear, ‘Computational Methods in Ordinary Differential Equations”, Prentice-Hall, Englewood Cliffs, NJ, 1971. (24) J. M. Haile and S. Gupta, J . Chem. Phys., 79, 3067 (1983). (25) R. D. Swindoll and J. M. Haile, J . Comput. Phys., 53, 289 (1984).
I
I
1 I
I -
Figure 1. Average singlet distribution functions pi(?) for each methyl and methylene group on a chain. Bars indicate estimates of statistical uncertainty, based upon scatter in the data at the locations of the bars.
pi(r) are shown in Figure 1. Bars on the curves in the figure indicate our estimates of the statistical uncertainties as reflected in scatter in the data. We make the following observations from Figure 1. (1) Only tails (i = 13) and CHz groups with i = 2 have measurable probabilities for their centers being found at or beyond the head-group sphere (ri 1 16 A). Thus, a chain’s central groups are nearly always within the interior, even when the tail is near the outer shell. However, only a small fraction of molecules have end portions of their chains ( i > 6) farther from the micelle center than head portions (i < 6). This tends to substantiate that aspect of the Dill-Flory m ~ d e lwhich ~ , ~ assumes that as i increases, groups either remain in a particular shell of several angstroms thickness or are directed away from the head toward inner shells. (2) The peaks in the p i ( r ) broaden as i increases, Le., as one moves from head to tail. Moreover, the peak heights decrease from i = 1 to 8 and then increase from i = 8 to 13. The decrease from i = 1 to 8 is probably because of the fixed head groups. The groups near the tail have greater asymmetry in p,(r) primarily because the presence of gauche bonds at small i has a strong effect on the radial positions of groups near the tail. The closer a group i is to the tail, the more possible combinations of gauche and trans bonds between groups 2 and i, hence the more radial positions accessible to group i. Preliminary simulation results for a second micelle model show that the p I distributions for groups i = 2-8 are broadened and lowered when the constraint of fixed heads is relaxed, however, no large changes are found for the distributions of the other groups. (3) The first nonzero values on the small-r side of the p,(r) curves are separated by 1.2-1.3 A, which compares with separations of 1.28 A for all-trans molecules lying along a micelle radius to a head group. The sharpness in the rise from zero of the pi for i = 2-4 is due to the head groups being fixed, but their spacing is not. In contrast, the peak heights are at different separations, bbth because gauche bonds change the radial positions of groups and because vectors along the chains are generally not aligned along the head-group radius vector. (4) No center of a tail group was found closer than 0.8 A to the micelle center. This does not imply a void at the center, since the effective methyl radius is more than twice 0.8 A. While a low probability of tails at the micellar center should be an excluded volume effect of three or four tail groups competing for the same location, the fact that none succeed in attaining the center appears to be a result of the radius chosen for the head-group shell. The point in item (4) is illustrated in Figure 2 which shows the interiors of the hemispheres produced upon slicing the micelle in two. The figure is at an arbitrarily chosen instant in time, 106 ps into the equilibrium segment of the simulation. Figure 2 confirms that little free volume exists at the micelle center (the groups are drawn with radii = u), that tail groups are concentrated at the center, and that tails may also be found at radial positions away from the center. The difference in tail concentration in the two halves of Figure 2 is an instantaneous effect; the time average, tail-group, center of mass coincides with the micelle center.
Intramicellar Structure via Computer Simulation
The Journal of Physical Chemistry, Vol. 88, N o . 25. 1984 6365
in the Lennard.Jones (6.9) potential
for each group
The results from these inte
average of the 13 singlet functions p,(r) gives the overall local density p(r) across the micelle This local density increases monotonically from zero at the micelle center to p(r) ‘v 0 0028 chain/,&’ at r ‘v 7 A, is essentially constant in the range 7 < r < 9 A, and then decreases in an oscillatory fashion from r IV 9 A to the head-group sphere at r = 16 A. These ascillations occur at radii associated with groups 2-7 and are thus probably caused by the head groups being immobile Figure 3 provides another description of the radial distribution of CH, and CH, groups The figure shows the conformations of eight molecules projected onto the ~4 plane, where r and 4 are two of the usual spherical coordinates (in a spaced-fixed frame): r, is the radial position of group I relative to the micelle center, and 4) is the group’s azimuthal angle. The eight molecules are those whose heads are fixed on an equator of the model micelle. The projections are shown at 106 ps into the equilibrium run and subsequently at intervals of 2 ps Three immediate observations can be made from a) The chains are clearly not in all-trans conformati a few conformations occur in which the chains ’curl hack” toward
mation at equilibrium. The average number of gauche bonds 16 therefore three per molecule. Whether this would he changed by a more realistic model IS uncertain. In the Weber and Ryckaert-Bellemans simulations of liquid alkanes, the gauche fraction is more than 40%, compared to the ideal-gas value of 34% In contrast, the bilayer simulation of van der Plwg and Berendsen gave 25%gauche bonds. The effect of density is likely to be small, according to both Weber and Ryckaert-Bellemans, so the straightening of chains in surfactant aggregates is probably real Table I1 presents the probability of finding a particular number of gauche bonds on one chain, determined by sampling our entire equilibrium simulation. The most probable number or gauche bonds per molecule is three, which agrees with tbe average number Only 2%of the molecules were found in the all-trans conformation, 75% had two, three, or four gauche bonds, and no molecules were found hanng eight or more gauche bunds Further, we computed the probability for any of the ten bonds on a chain to be gauche Within our statistics, all bonds were equally likely to he gauche van der Ploeg and Berendsen also found a uniform distribution of gauche bonds along a chain in their bilayer simulations The radial dependence of the gauche fraction was also determined, No direct one-to-one correspondence exits between the
6366 The Journal of Physical Chemistry, Vol. 88, No. 25, 1984 TABLE III: Probability of Finding Two Gauche Bonds on the Same Chain, with the Two Gauche Bonds Separated by Different Numbers of Bonds
separation in bonds 0 1
2 3 4
probabilitv
separation in bonds
probabilitv
0.147 0.202 0.156 0.148 0.117
5 6 7 8
0.094 0.069 0.047 0.020
’- SHELL I
OF FIXED HEADS’ R=168
.\
,\
/
/ \
Figure 4. Orientations of chains relative to the micelle radius through the head group. On the left is the average orientation of all molecules found in the all-trans conformation. On the right is a “typical”molecule
that, at one instant in time, has several conformational features similar to those of an *average”molecule. The group diameters used in this figure are 0.670 (2.34 A) to more clearly display the molecular conforma tions. position of a group or bond on a chain and the radial position of that group or bond in the micelle (except, of course, for the fixed head groups of our model). For bonds in the range 6 < r < 14 A, the gauche fraction was essentially constant at the average value of 29%. For both r < 6 A and r > 14 A, the gauche fraction increases to loo%, but this should not be. taken too seriously since the statistics for sampling dihedral angles are poor at these extreme values of r. Table I11 presents a pair probability function for gauche bonds on a chain; Le., given that a molecule has two or more gauche bonds, what is the probability of a particular number of bonds (trans or gauche) separating two gauche bonds on the chain? The table indicates the most probable separation is one bond between two gauche bonds, although the probabilities for separations by zero, two, or three bonds are also high. Thus, 65% of all gauche pairs on a chain are separated by no more than three intervening bonds. The third basic question we considered was the average chain conformation relative to the micelle radius vector through the head group. While only 2% of the molecules were found to be all-trans,
Haile and O’Connell even these are not necessarily aligned along a micelle radius. As shown in Figure 4, the average orientation of the all-trans molecules is such that the tail group is about 8.4 A from the micelle center. Also shown in Figure 4 is the orientation of a “typical” molecule. The molecular conformation shown in the figure is one actually observed at one instant of fime in the simulation. This typical molecule has several features characteristic of the “average” molecule: (a) it has three gauche bonds; (b) the distance from the tail to the micelle center is 9.3 A (cf. R I 3= 9.27 A in Table I); (c) the end-to-end distance of the chain is 13 A (cf. a most probable end-to-end distance of 12.8 A (though the average is 11.6 A)); (d) the angle subtended by the radius vector to group 13 and the radius vector to the head is 54O, while the average angle for all molecules is 42’.
Conclusions We believe that although the model micelle studied here is simplified, many of its qualitative aspects are realistic. (1) Unlike a traditional description of intramicellar structure, we find only a very small fraction of chains are in an all-trans conformation, and even these are not, on the average, aligned along a micelle radius. While a concentration of tails exists near the micelle center, the density of groups is not higher in the center thaii in outer shells. (2) A significant fraction of bonds are in gauche conformations. At the density studied here, the average gauche fraction of the model is 29%. Because this was an initial simulation, we purposely chose to use a smaller number of chains than occurs in a dodecyl surfactant micelle, so that a tolerable amount of computer time would be needed for the simulation. The effect of higher density is not yet known. (3) N o simple correlation exists between position on a chain and radial position in a micelle. Simulations now in progress do not have the obvious limitations of this first simulation, namely the fixed head groups and the low chain density. Other improvements could include more realistic chain-solvent interactions and better inter- and intramolecular potentials models for the chain-chain interactions. Acknowledgment. This simulation was performed on the IBM 370/3081 at Clemson University using computer time made available through the College of Engineering. Graphical analyses were done on the University’s DEC VAX 11/780 using a hidden-line program originally written by Nelson Max at Lawrence Livermore Laboratory, with modifications and additions by Dr. S . M. Thompson, School of Chemical Engineering, Cornel1 University, and by S.Gupta at Clemson. Partial financial support for J.M.H., including a Presidential Young Investigator Award, is provided through the Chemical and Process Engineering Division of the National Science Foundation.
