Langmuir 1993,9, 916-926
916
Molecular Dynamics Simulation of an Aqueous Sodium Octanoate Micelle Using Polarizable Surfactant Molecules John C. Shelley,+Michiel Sprik,* and Michael L. Klein'lg Department of Chemistry, University of British Columbia, Vancouver, British Columbia V6T lZ1, Canada, Zurich Research Laboratory, IBM Research Division, CH 8803, Riischlikon, Switzerland, and Department of Chemistry, University of Pennsylvania, Philadelphia, Penmylvania 19104-6323 Received June 1, 1992. I n Final Form: December 23, 1992 The inclusion of dynamicinduced polarization in simulationsof aqueous ionic systems can have important consequences. As a step toward performing analogous studies on aqueous surfactant systems, we have completed initial studies of the system sodium octanoate/water in which only the surfactant anions are polarizable,polarization effects in the solvent being treated via a mean-field model. Two simulationswere performed: one of an aggregate of 15 sodium octanoate monomers and the other of a single octanoate ion in solution. The present simulations are compared to earlier results for nonpolarizable models. The behavior of the monomers is not appreciablychanged by the addition of polarization. However, general structural characteristics of the micellar aggregate such as volume, shape, and hydration are different. Significantly,with the current parameterization,no contact ion pairs of the type found in the nonpolarizable model were present when the anion headgroups were allowed to polarize,a result that accords better with the interpretation of experiment. 1. Introduction
Computer simulations of atomic and molecular systems using Monte Carlo or molecular dynamics methods that employ classical Hamiltonians are widely used. The potential energy and force calculations usually involve using a completely predetermined potential energy function which in most cases relies on the assumption of pairwise additivity; i.e., explicit many-body effects such as the polarization (induction) interactions are not included. For macromolecules, the pairwise additive model of the potential is augmented using intramolecular bond, bond angle, and dihedral angle interactions, among others. The great advantage of this approach is that it is relatively inexpensive and in many cases the only practical method currently availablefor simulatinglarge systems. However, some many-body effects may not be modeled accurately enough using effective pair potentials.lP2 For example, in aqueous ionic systems,the major many-body contribution is that involving induced dipole interactions. A number of schemes have been proposed to handle the induced polarization of ~ a t e r . ~ -We l ~ will focus on the recently +
University of Birtish Columbia.
$
University of Pennsylvania.
* IBM Research Division.
(1) Allen, M. P.; Tildesley, D. Computer Simulations of Liquids, Clarendon Press: Oxford, 1987. (2) Price, S. L.; Stone, A. J. J . Chem. Soc., Faraday Trans. 1992,88, 1755. (3) Barnes, P.; Finney, J. L.; Nicholas, J. D.; Quinn, J. E. Nature 1979, 282. 459. ~
~~
(4jCaillo1, J. M.; Levesque, D.; Weis, J. J.; Kusalik, P. G.; Patey, G. N. Mol. Phys. 1985,55, 65. (5) Caldwell, J.; Dang, L. X.;Kollman, P. A. J. Am. Chem. SOC.1990, 112,9144. (6) Rullmann, J. A. C.; van Duijnen, P. Th. Mol. Phys. 1988,63, 451. (7) Sprik, M.; Klein, M. L. J . Chem. Phys. 1988,89, 7556. (8) Wallqviest, A.; Ahlstrbm, P.; Karlstrbm, G. J . Phys. Chem. 1989,
-.w. -,1Rd9. -(9) Ahlstrbm, P.; Wallqviest, A.; EngstrBm, S.; Jonsson, B. Mol. Phys. 1989,68, 563. (10) Sprik, M.; Watanabe, K.; Klein, M. L. J . Phys. Chem. 1990, 94, RAR.? _"". (11) Niesar, U.; Corongiu,G.;Clementi,E.;KneUer,G.R.; Bhattacharya, D. K. J . Phys. Chem. 1990,94,7949. (12) Cieplak, P.; Kollman, P.; Lybrand, T. J. Chem. Phys. 1990, 92, 6755. (13) Straatsma, T. P.; McCammon, J. A. Mol. Simul. 1990,5, 181. ~
developed extended Lagrangian scheme7J5J6that employs additional classical degrees of freedom to handle the induced polarization. With suitable parameterization one can reproduce the structural and dynamic properties of water.16 The scheme has been applied to simple ions in ~olution.~~J7J8 The present paper describes an attempt to extend this approach to more complex polyatomic ions such as those in ionic surfactant system^.^^-^^ We focus particularly on micelles of sodium o c t a n ~ a t e . ~This ~*~~*~~ system is not only of intrinsic interest since experience gained in this current effort may be useful in developing polarizable models for more complex biological systems including phospholipids, peptides, and proteins. Questions have arisen about the validity of using pairwise potentials to describe crucialaspecta of surfactant systems. The headgroup of a surfactant anion is quite polarizable and, as a result, may be significantly influenced by the presence of nonsymmetric charge distributions in ita environment, whether there is a counterion nearby or not. Also, in most simulations heretofore, the effectivedielectric constant used to treat the hydrocarbon micelle core is 1 as opposed to the experimental value of roughly 2. Since chargeinteractions between headgroups acrossthe micellar aggregate may be important for determining structural stability, proper treatment of the dielectric nature of the hydrocarbon core may be important. We have chosen to carry out an initial study of the sodium octanoate system including both intramolecular and intermolecular induced polarization in the amphiphiles to investigate the changes that occur. Since the (14) Motakabbir, K. A.; Berkowitz, M. L. Chem. Phys. Lett. 1991,176, 61. (15) Sprik, M. J . Phys. Chem. 1991, 95, 2283. (16) Sprik, M. J. Chem. Phys. 1991,95, 6762. (17) Sprik,M. MetalsinSolution. In JournaldePhysiqueIV, Colloque C5; Damay, P., Leclercq, F., Eds.; 1991; Vol. 1, pp 99-102. (18) Karim, 0. A. Chem. Phys. Lett., submitted for publication. (19) Watanabe, K.;Ferrario, M.; Klein, M. L. J . Phys. Chem. 1988,92, 819. (20) Watanabe, K.; Klein, M. L. J. Phys. Chem. 1989, 93, 6897. (21) Wahnabe, K.; Klein, M. L. J . Phys. Chem. 1991,95, 4158. (22) Shelley, J. C.; Watanabe, K.; Klein, M. L. Int. J . Quantum Chem. 1990, 17, 103. (23) Shelley, J. C.; Watanabe, K.; Klein, M. L. Electrochim. Acta 1991, 36, 1729.
0743-7463/93/2409-0916$04.00/0Q 1993 American Chemical Society
Aqueous Sodium Octanoate Micelle
Langmuir, Vol. 9, No. 4, 1993 917
number of polarization sites is small compared to the total number of interaction sites in the system, the cost in computer time is much less than if one introduced polarization sites for the entire system. In the near future we plan to include the polarizability of the solvent. The octanoate system is also appropriate for this initial study because on one hand it is small enough to permit lengthy calculationsand on the other hand it is well characterized by our earlier work.1g-21.23
was employed for the molecular dynamics simulations at a temperature of 298 K. For the simulation of the octanoate monomer in water we used DzOinstead of HzO which permitted a larger time step (3.5 fs as opposed to 2.5 fs). The POL model is based on a scheme developed previously to study water and ionic solutions using an extended Lagrangian ~cheme.~J”l~ In this approach, additional charge sites are added to the moleculesand the values of the charges assigned to these sites are treated as 2. Molecular Dynamics Calculations classical dynamic variables which evolve (like any other degree of freedom) in response to the local electric field. The simulation method is essentially the same as that For a detailed description for the case of water the used previously.mp23However, in the current study induced interested reader should consult ref 15. For polarizable polarization of the surfactant ions necessitated some ions, the model is similar but not identical. A preliminary corresponding changes to the potential parameters of report of an application to the aqueous C1- ion can be surfactant counterion interactions. We begin by briefly found in ref 17. In the present study, we will limit the reviewing the original nonpolarizablemodel (NP)and then discussion to a brief overview. The induced polarization describe how it was changed to arrive at the current of monoatomic ions such as halides can be handled using polarizable surfactant model (POL). two additional interaction sites. The pair of extra sites, The NP model is based on pairwise additive site-site plus the nucleus in between, form a rigid triatomic nonbonded interactions consisting of Lennard-Jones pomolecule. The middle site carries the total ionic charge tentials along with Coulombic interactions between point and is also the location where the short-range (Lennardcharges on some of the atoms. The SPC/E effective pair Jones) interaction is centered. The charges on the outer potential24 was used to describe the water molecules. The sites are equal in magnitude but of opposite sign and vary potential parameters for Na+,which carried a full l e charge, with the applied electric field in a dynamic manner. As were taken from the literature.25 Hydrogen atoms were a result, the ion acquires an instantaneous dipole moment not explicitly represented in the hydrocarbon tails of the with a magnitude determined by the value of the dynamsurfactant. Instead, each CHZ and CH3 group was ically fluctuating charges. A change of orientation of the described by one mass and force center referred to as a dipole is brought about by the rigid “triatom” pivoting pseudoatom. The carbon atoms or pseudoatoms have been about its center. The finite-size linear molecular ion numbered 1-8 starting from the headgroup. The latter behaves effectively as a polarizable monoatomic ion, consisted of two explicit oxygen atoms and the C1 atom provided the bond length is small in comparison to the which carried charges of -0.7e and 0.4e, respectively, for dimensions of the ion. a net charge of -le. Each of the headgroup atoms also The value of the polarization at a given time is carried a Lennard-Jonesforce center. Lorentz-BerthelotZ6 determined from regular second-orderequations of motion combining rules were used to determine the nonbonded derived from the extended Lagrangian. The principle interactions between force centers of different types. The underlying extended Lagrangian schemes is that the Lennard-Jones potential parameters for the octanoate ion fictitious dynamics of the pseudo degrees of freedom (the were taken from ref 27. dipoles) approaches the adiabatic limit of the instantaThe intersite distances of the SPC/E water molecule neous reponseto the slowlychangingatomicconfiguration. were fixed using the method of Bond lengths This imposes certain restrictions on the specification of in the octanoate anion, including the geometry of the the mass of the two end sites and their charge. These headgroup, were fixed using Shake.29 Intramolecular parameters control the response time of the polarization b~nd-bending,~~ and nonbonded potentials and should be chosen sufficiently small to ensure sepa(between sites further away from each other than three ration of time scales between the rapid fluctuations of the bonds)were used. There was no torsionalpotential applied dipoles and the true physical motion of the system. A to the rotation of the headgroup about the c1-C~bond. condition of adiabatic locking is established if the polarA 12-A cutoff distance was used for the nonbonded ization fluctuationsare also cooled,Le., if their temperature interactions in the micelle simulation, whereas for the is constrained to a value low compared to the thermosingle monomer in water, the cutoff was 10 A. In the dynamic temperature of the solvent and solutes. In ref micelle simulation, triclinic periodic boundary conditions 15 it was shown that static properties such as energetics were used with the three box vectors 32.04 A long and the and structure are not affected too much if the adiabatic three angles all 109.4’7’. The use of such a cell reduces the condition is not satisfied and the system is in thermal number of water molecules needed to keep the periodic equilibrium. However, it was found that cooling and images of the micelles apart. An Ewald method’ was used tuning of the dynamics of the polarization are absolutely to calculate the long-range electrostatic contributions to in order to obtain the correct dynamical the potential energy. The canonical (N, V,T,e n ~ e m b l e ~ l s ~ necessary ~ behavior. In the present simulation, we employ the nonequilibrium adiabatic scheme and keep the temper(24) Berendsen, H. J. C.; Grigera, J. R.; Straatsma,T. P. J.Phys. Chem. ature of the polarization at 5 K using NosbHoover 1987,92,6269. (25) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. therm0stats3~for both the size and orientation of the Chem. Soc. 1984,106,903. dipoles. Because of the strong coupling of the induced (26) Hirshfelder, J. 0.; Curtis, C. F.; Bird, R. B. Molecular Theory of momenta to the librational motion of the water molecules, Gases and Liquids; John Wiley: New York, 1954. (27) Jorgensen, W. L. J . Am. Chem. SOC.1981, 103, 335. a complementaryset of thermostats is applied to maintain (28) Ciccotti, G.; Ferrario, M.; Ryckaert, J. P. Mol. Phys. 1982, 47, the temperature of their configurationaldegreesof freedom 1253. at the appropriate value. (29) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput.Phys. 1977, 23, 327. Another important issue in the modeling of polarization (30) van der Ploeg, P.; Berendsen, H. J. C. J. Chem. Phys. 1982, 76, concerns the damping at short range. The experience with 3271. (31) Nose, S . Mol. Phys. 1984,52, 255. (32) Nod, S. J . Chem. Phys. 1984,81, 511.
(33) Hoover, W. G. Phys. Reo. A 1985, 31, 1695.
918 Langmuir, Vol. 9, No. 4, 1993
water has convinced us that some form of control of the electrostatic induction by polar molecules must be included, particularly if a pair is joined in a direct hydrogen bond. In our representation of polarization, short-range damping can be readily implemented by replacingthe point charges with Gaussian charge distributions at the polarization ~ i t e s . ~ JThe ~ Jwidth ~ of these Gaussiansis adjusted so that known properties of the system,such as the average net dipole moment of the molecules in question, are reproduced. The need for a simple and stable short-range damping is, in fact, the main motivation for modeling the induction in a monoatomic particle by a three-site rigid rotor with fluctuating charges. As an alternative one could consider a Drude oscillator consisting of two particles of fixed charge connected by a spring. However, in the rotor model the damping is strictly independent of the magnitude of the dipole moment, while in the spring model the variable separation between the positive and negative charges leads to occasional but potentially catastrophic instabilities. In the POL model, the basic structure of the octanoate anion was not changed except that each atom or pseudoatom was provided with a polarizable center of the type described above. The distance between the variable charges at the end points of the rigid rotor was set at 0.30 A. The middle site which carries the permanent charge was fixed at the originalmass center. Each of the dynamic charge sites was assigned the mass of a H atom, and the mass of the original atom or pseudoatom was correspondingly reduced by 2 amu. The generalized mass used for the oscillating charge degree of freedom was 0.5 m u . Our choice of inertial parameters increases the time scale for the dynamics of the induction to close to the librational time scales of water. This is a compromise between a full separation of time scales required for a rigorous adiabatic coupling of the polarization to the configurational degrees of freedom and a value of the molecular dynamics time step, which enabled us to follow the evolution of the polarizable system over periods of several hundred picoseconds. The time step of 2.5 fs that was used is, in fact, typical for the simulation of aqueous systems. Since the time step was limited by the polarization, we did not run the simulation with deuterated water. The carboxylate headgroup is the most polar part of the anion, and the treatment of the intraionic coupling of the permanent charges on the C1,01, and 02to the inductions is the next critical step in the modeling of the polarization of the octanoate ion. All interactions, charge-dipole and dipole-dipole, of the C1 carbon with the two oxygens are excluded in our scheme while the mutual interactions of the oxygen atoms are retained. This distinction is based on the following argument. The C1-0 bond is strongly polar, and we assume that polarization effects across such a bond have already been taken into account in the choice of assigned permanent charges. The oxygen atoms, on the other hand, are not linked by a direct chemical bond and hence are treated similar to any other nonbonded pair. For the carbon chain the situation is somewhat ambiguous. Sincethese bonds are homopolar and covalent, it was decided to permit induction coupling, including the charge-dipole interaction across the bond. The justification for not inhibiting the effect of the C1 charge on C2 is that this charge is part of the representation of the headgroupand is, therefore,unrelated to the essentially homopolar C& bond. Once the pattern of allowed intramolecular interactions is defined, the model parameters can be specified. The permanent charges on the 0 and C1 sites are the same as in the nonpolarizablemodel. The site polarizabilities, ai, were determined via a decompositionof the polarizabilities
Shelley et al.
Figure 1. A gas-phase polarizable (POL) octanoate anion at very low temperatures depicted along with pairs of induction sites for each atom or pseudoatom. Positive (light) and negative (dark) induction sites appear as spheres with radii proportional to the size of the induced charge. For reference, the induced charges on the oxygens have a magnitude of about 0.6e.
of simple alkanes and common functional groups.26 This gave polarizabilities of 1.16, 1.87, and 2.24 A3 for the C1, methylene (CH3, and methyl (CHd groups, respectively. We chose to use 1.5 A3 for each of the headgroup oxygen atoms. The resulting distribution of self-induced dipole moments of an isolated surfactant ion is shown in Figure 1. The dominant feature that characterizes the model and that will alsohave implicationsin solutionis the sizable dipoleon an oxygen atom ( 1D) arisingfrom the coupling to the permanent charge on the other headgroup oxygen. If the equivalent model is applied to a neutral alcohol or carboxylic acid molecule, it is found that the total induced dipole moment of the molecule points in a direction approximately opposite to that of the dipole moment associatedwith the permaent charges,i.e., the "bare" dipole moment of the nonpolarizable model. The permanent dipoles in the class of potentials used here are, in general, larger than the experimental gas-phase values. This observation suggested a convenient procedure for fixing the value of the remaining model parameter, namely, 4, the width of the Gaussian used for the induction sites. By increasing this length parameter until the experimental gas-phase values were recovered, we obtained a value of 4 = 0.7 A. This value is roughly consistent with that found to be appropriate for polarizablewater in previous work.l6 Values for the dipolesinduced on each site in the octanoate ion are listed in Table I11 based on this parameterization. The final stage of parameterizing a polarizable model focuses on the interaction of the sodium and octanoate ions in the gas phase. For initial trials, we used the NP potentials without change along with the induction model describedabove. In accordancewith previousexperience,17 such a "unrenormalized" potential yields a Na+ binding energyfar too large. The reason is that to some extent the original Na+ Lennard-Jones parameters were chosen to effectively include the charge-induced dipole interaction that the ion induces in neighboring species. In the POL model such interactions are explicitly accounted for. Hence, we were double-counting this contribution to the potential energy. An obvious way to correct this situation is scaling back the permanent charges. However, as indicated above, in the present study we opted for retaining the NP charges and have modified only the Lennard-Jones parameters. Adoptingthis constraint was largely a matter of convenience. It is certainly one of the elements in our parameterizationthat has to be reexamined,in future work. A more systematic parameterization scheme, perhaps related to that developed by St0ne,3~.~~ is clearly needed! In a manner consistent with the Lorentz-Berthelot26 combining rules, the Lennard-Jones c and e for all of the interactions of atoms or pseudoatoms in the octanoate with the Na+ were collectively adjusted to get the same geometry and roughly the same binding energy for the Na+ to the octanoate headgroup as were obtained for the
-
(34) Stone, A. J. Mol. Phys. 1985,31, 1695. (35)Stone, A. J.; Alderton, M. Mol. Phys. 1985,56, 1047.
Aqueolls Sodium Octanoate Micelle
Langmuir, Vol. 9, No. 4, 1993 919
Table I. Lennard-Jonee Parameters for Nat.CH2(CH2)&02- Interactions in the Polarizable Model
0 CI
3.107 3.184
21.18 10.762
CHz CH3
3.557 3.496
22.458 28.28
NP model. The electrostatic energy of the POL model is offset by a positive polarization energy (p2/2cr). As noted earlier, the octanoate ions internally induce dipoles with significant polarization energies, even in the gas phase. The difference in polarization energy of sodium octanoate in a bound and dissociated state is 13 kJ/mol and may be much larger for other accessible configurations. This is another illustration of how the energy balance can be dramatically changed by polarization. Actually, the Na+ ion is still bound more strongly by about 12 kJ/mol in the POL model, but further modifications to the LennardJones parameters seemed unreasonable. The final parameters are listed in Table I. The modifications are consistent with those employed for the C1- ion interacting with polarizable water.1° We did not adjust the corresponding potential parameters for the water-octanoate ion potentials. The resulting increase in the wateroctanoateinteractionsbrought them into better agreement with quantum mechanical calculations and an earlier potential model (see ref 36 and references therein) for the water-octanoate geometries found in previous simulat i o n ~as~well ~ *as~the ~ current one. Differences between the associationof water or ions with the headgroup for the NP and the POL models should not be taken as being solely due to the inclusion of polarization effects. The effect of the renormalization of the atom-atom potential is at least of equal importance. Initially,the complexityof the code restricted us to using a cubic system, so a well-equilibrated configuration obtained from our earlier simulation23 with triclinic periodic boundary conditions could not be used directly. Instead, we periodically replicated such a configuration and cut out a cubic box containing the 15 surfactant ions from one micelle only and randomly removed any excess Na+ ions. The resulting system contained 1204 water moleculesas compared to the original system of 718 water molecules. During annealing, the volume was adjusted by a small amount to get the expected volume, 40 000 A3 for the composition of this system (the length of the box vectors was 34.2 A). The resulting concentration of surfactant was 0.62 M. The monomer simulation was a continuation of the earlier and consisted of 1octanoate ion in 384 water molecules in a triclinic box with all box vectors 28.5 A and all angles 109.47'. T w o simulations of the micelle system were initiated: one for the POL model and another using the NP model, for reference. Each system was propagated for an initialization period of roughly 30 ps. During this time essentially all of the contact ion pairs in the POL system broke up while there was little change in the NP system. Then, the POL system was run for a total of 450 ps on an IBM 3090 6005 at the Cornell Supercomputing Center (taking full advantage of the parallel and vector facilities there). The monomer system was run for a total of 1025 ps on our local Convex C210. 3. Micelle Results
Some representative system properties for the POL model, calculated at regular intervals during the simu(36)Jorgensen, W. L.; Gao, J. J. Phys. Chem. 1986,90, 2174.
lation, are given in Table 11. While considerable fluctuations are present, most properties do not appear to exhibit significant long-term trends, which suggests the system has settled down. The exception to this is possibly the ion association data. Further analysis, described later in this paper, shows no long time changes in ion association are found after 150 ps. All averages and time correlation functions reported in this paper therefore span the remaining 300 ps of the simulation. Numerical and visual checks on a graphics workstation confirm that the micelle stays together. Also no monomers left the micelle during the run. Configurations near the beginning and end of the 300-ps trajectory are depicted in Figure 2. In view of the obvious similarities of the two parts of the figure, it is meaningful to speak of the micelle as an "apparently stable aggregate" and to calculate properties based upon this concept. For the purposes of the presentation and analysis of the results we now consider,in turn, the waterhydrocarbon interface, the micelle shape and size, Na+ ion association, and the properties of the octanoate ions. 3.1. Water-Hydrocarbon Interface. Figure 3 gives the average radial densities of carbon atoms, headgroup oxygen atoms, Na+ counterions, and water oxygens with respect to the center of mass of the micelle. The POL simulation results are shown (bold lines) along with those for the nonpolarizableNP micelle.23The volume sampled for distances close to the center of mass of the micelle is small. Hence, the uncertainties in the hydrocarbon densities are large,and the differencesbetween the results for the two simulations are likely not significant in this region. Water penetration into the micelle is enhanced in the simulation with polarization; however, there is still a region about 5 A in radius that is essentially devoid of water. There are on average 93 water moleculesin contact (within 4.0 A) with the hydrocarbon tails as compared to 79 for the nonpolarizable system. Since more than half of these water molecules have neither a headgroup nor Na+ as a nearest neighbor, there is a water-hydrocarbon interface. It would be instructive to characterize such an interface on the molecular level. Since the data in Figure 3 are averaged over all times and angles, the breadth of the region of overlap of the hydrocarbon density and the water density may be in part a result of the nonspherical shape of the micelle (recall Figure 2) in addition to fluctuations about the average volume and shape. To try to learn more about the nature of the interface, it is perhaps appropriate to attempt to define a "local" surface which represents the point at which the local hydrocarbon density achieves a specified value. Straightforwardly calculating the local densities from the instantaneous positions of the atomic nuclei and using them to define and characterize a local surface would be unwieldy and likely yield noisy results. Instead, an effective local density may be calculated by placing a broad, normalized Gaussian on the nucleus of every carbon atom ((2248) in the hydrocarbontails. The sum of all such Gauasiam yields a function which can be used as a measure of the local hydrocarbon density. The use of Gaussians effectively allows for some of the thermal motion of the atoms and yields a smoothly varying measure of the local density. A contour of this function at some test value may then be used to define a reference surface. Average densities of the various components of the system may be calculated relative to this reference surface, yielding information about the location of components of the system relative to the interface. A series of tests were performed to determine the sensitivity and dependence of the estimated breadth of the interfacial region to the width of the Gaussians and the value of the density contour employed. The breadth
920 Langmuir, Vol. 9, No. 4, 1993
Shelley et al.
Table 11. Results for the POL Model of Sodium Octanoate ~
Na+ region tlps
0 50 100 150 200 250 300 350 400 450
(PE)"/(kJmol-')
R(Cl)*/A
R(Cd*/A
percent gauche'
Id
le
fred
water contact.@
-57.57 -57.56 -57.56 -57.57 -57.58 -57.60 -57.60 -57.60 -57.60
10.5 10.9 10.2 9.8 10.3 9.8 10.9 11.0 10.7 10.2
7.1 6.0 7.0 7.1 7.3 7.3 8.0 7.5 6.6 7.9
15 23 25 24 24 23 21 21 12 29
0 1 1 0 0 0 0 1 0 0
7 4 2 7 8 6 8 8 7 7
8 10 12 8 7 9 7 6 8 8
90 85 104 93 83 82 102 100 101 98
*
Averaged over the previous 50 ps. Instantaneous average over all chains. Instantaneous average over all chains and dihedral angles. Ions within 4.0Aof a C1atom. e Ions between 4.0 and 6.55A of a C1 atom. f Ions further than 6.55A from a C1 atom. 8 Water molecules within 4.0A of any methylene or methyl pseudoatom a t that instant. d
of the interfacial region was arbitrarily taken to be the distance between points at which the hydrocarbon and water densities cross through the value 0.003 A-3 (about 10% of their bulk densities). The estimated width increased monotonically with the width of the Gaussians. However,for Gaussians of width from 1to 2 times the size of the water-CHz Lennard-Jones u parameter (about 3.67.2 A), the estimated interfacial width only changed from roughly 2 to 3 A. Such a width effectively determines a surfaceby averagingover the locationsof several molecules. The major effect of changing the densities used for the test surface within the range of 10-40% of the bulk hydrocarbon density was to shift all of the average density profiles of the hydrocarbons, water molecules, and headgroup oxygens by roughly the same distance relative to the reference surface without significantly affecting the width of the interface. Figure 4 shows roughly half of the system at the end of the simulation along with the corresponding reference surface for Gaussians of width 1 . 5 and ~ a test value at 25% of the bulk hydrocarbon density. With the same parameters, the average densities of the tail carbons (C2-C*), headgroup atoms (C1 and the two oxygen atoms), Na+ counterions, and water oxygen atoms were determined as a function of the distance from the local surface by averaging over the simulation (see Figure 5 ) . Here, the interfacial region is delimited by vertical lines and is roughly 3 A wide. Essentially, all of the head roup atoms lie inside this interfacial region or within 2 of it in the aqueous phase, and most of the Na+ ions lie in the aqueous phase, The large narrow peaks in the Na+region are due to the large statistical noise present in this method. The narrowness of the interfacial region and the placement of the headgroups just outside this region are in good agreement with the interpretation of experiments for similar systems.37~38 The present results are consistent with those obtained using similar methods on SDS.22 The estimated width of the interface is similar to, but somewhat narrower than, results obtained in other theoretical studies.19*20~39~40 However, since those estimates were not based on a local surface, some diferences are to be expected. 3.2. Micelle Size and Shape. The interface between the hydrocarbon part ofthe micelleand the water is narrow. Accordingly, it is appropriate to talk about the micelle's size and shape. The approach used in previous studies20 was to define a micellar volume
x
s2 = (4/3)dRHG3)
(1)
(37) Cabane, B. J. Phys. 1981, 42,847. (38)Stigter, D.J. Phys. Chem. 1974, 78, 2480. (39)Koopal, K.K.;BBhmer, M. R.; Lyklema, J. J. Phys. Chem. 1991, 95, 9569.
(40)Jonsson, B.; Edhom, 0.;Teleman, 0. J . Chem. Phys. 1986, 85, 2259.
using the headgroup positions at each time step. Figure 6 showsthe evolution of s2 over the course of the simulation determined by this procedure. The average value of tl for the POL model is 5600 A3 as compared to that for the nonpolarizable system, 5000 A3. If one calculates an average headgroup radius on the basis of the volume, it is 11.0 A as compared to 10.6 A for the nonpolarizable simulation. As in previous work,20we have also employed the instantaneous principle moments of inertia of the micelle, 11,12, and 1 3 , to monitor its shape. The evolution of these moments during the simulation is depicted in Figure 7. As before there are significant fluctuations, but two of the moments of inertia are consistently much larger than the third, indicating that the micelle maintains a prolate shape. The average ratios of these two large moments to the smaller one are 1.94 and 2.10 versus 1.19 and 1.37 for the nonpolarizable system, indicating that the polarizable micelle is more distorted from a spherical shape. The time autocorrelationfunctionof the deviations from the average volume and the average ratios of the two largest moments of inertia is quite noisy. The time constants for the relaxation of fluctuations of the volume and ratios of inertia are estimated to be roughly 50 ps which is longer than those for the nonpolarizable ~ystem.~3 Becausethe polarizablemicelle shape deviates significantly from a sphere, the overlap in the hydrocarbon and water densities extends over a wide range (recall Figure 3). The nonspherical nature of the polarizable micelle is clearly evident from the shape of the reference surface shown in Figure 4 and the instantaneous configurations shown in Figure 2. 3.3. Counterion Association. In the nonpolarizable the radial distribution function of the Na+ with respect to the headgroup C1 atom was sufficiently structured (see the dashed curve in Figure 8) to classify the ions into three cate ories, R c ~ - N
