J. Phys. Chem. 1994, 98, 712-711
712
Molecular Dynamics Simulation Study of an n-Decyltrimethylammonium Chloride Micelle in Water+ Josef Biickers and Jiirgen Brickmann' Institut fZr Physikalische Chemie I , Technische Hochschule Darmstadt, Petersenstrasse 20, 0-64287 Darmstadt, Federal Republic of Germany
Philippe Bopp Institut fiir Physikalische Chemie, Rheinisch- WestjZlische Technische Hochschule Aachen, Templergraben 59, 0-52062 Aachen, Federal Republic of Germany Received: September 9, 1993'
A molecular dynamics (MD) simulation study of a n-decyltrimethylammonium chloride micelle in water is presented. The model system contains 30 amphiphile molecules, 30 chloride ions, and 2166 water molecules. W e use the same model potentials as in a recent MD study of a n-hexadecyltrimethylammonium chloride monolayer in water (Biicker, J.; Schlenkrich, M.; Bopp, P.; Brickmann, J. J. Phys. Chem. 1992,96,9915). The main objective of this work is to analyze the shape and the structure of the micelle. A slightly prolate ellipsoidal shape of the micelle emerges from a simulation of 275 ps. It is found that the interior is completely "dry". This result is in contrast to monolayers-aqueous solution interfaces (Biicker, J.; Schlenkrich, M.; Bopp, P.; Brickmann, J. J . Phys. Chem. 1992,96,99 15). A molecularly sharp interfacial region appears, in agreement with experimental studies. The charges on the surface of the micelle are only partially neutralized by chloride counterions in the first solvation shells of the head groups. The area per head group on the micelle surface and the orientations of the alkyl chains as obtained from the simulation are in accord with experimental results.
Introduction Micelles are small aggregates of amphiphile molecules in aqueous solution. There are many experimental studies on the thermodynamic and other phenomenological properties (see, for example, ref 2). In contrast to this fact there is not much known on the microscopic structure of micelles. Because of fluctuations in size and shape it is difficult to establish the detailed structure from experimental studies alone. Most information is obtained from small angle neutron scattering In these studies a simple model is used which assumes a certain size for a micelle and wherein its shape is determined by the ratios of the principal moments of inertia. Geometrical information is obtained from light scattering experimentsa8N M R measurements9J0give evidence about the conformation and orientation of the alkyl chains. The molecular dynamics (MD) simulation method is well suited to investigate many particle systems from a microscopic point of view and so fill the gap between theory and experiment. This has been recently demonstrated in a paper by us1wherein on the basis of extensive MD studies of a n-hexadecyltrimethylammonium chloride (C16TAC) monolayer in water' we predicted the structural results of a neutron reflectivity study of monolayers of n-tetradecyltrimethylammonium bromide (C14TAB) and n-octadecyltrimethylammoniumbromide (C18TAB).11 The calculated density profiles, perpendicular to the monolayer, of water, head groups, and chains are in excellent agreement with the ones from the neutron reflectivity study. So far, only a few authors have simulated a micelle with an adjacent water phase.12-1s Micelles without an explicit water phase were more often The comparison of these studies shows that the water phase has a significant influence on the interfacial structure Author to whom correspondence should be addressed. Part of Dr.-Ing. thesis of Jmef Bkker, TechnischeHochschuleDannstadt, D17, 1993. 4 Present address: Max-Planck-Institut for Chemie (Otto-Hahn-Institut), AG Physikalische Chemic, Joh. J.-Becher-Weg 27, Universtitscampus, D-55128 Mainz, FRG. *Abstract published in Advance ACS Absrracts, December 15, 1993. f
0022-3654/94/2098-07 12$04.50/0
and on the shape of the micelle, whereas the interior of the micelle is not strongly affected by the solvent as expected. The model system in the present study comprises 30 amphiphile molecules, 30 chloride ions, and 2166 water molecules. The number of amphiphile molecules was chosen by extrapolating the results of a study by Imae et a1.,8 because there are no suitable experimental data on the shape and the size of a CloTACmicelle. These authors found aggregation numbers N of 84 for C16TAC, 62 for C14TAC, and 44 for CI2TAC. An uncertainty of about * 5 molecules can be assigned to these numbers due to the influence of the concentration of the micelles. The critical micelle concentration, cmc, of CloTAC is reported as 0.05 mol/L23 and 0.0614065 mol/L.24 Weinclude that themodel system as much water as possible because most experiments are done at concentrations near the cmc. Measurements of aqueous sodium chloridesolutions of C14TA with different chloride concentrations have shown that a spherical micelle builds up at chloride concentrationslower than 2.7 mol/L.8 The chlorideconcentration in the simulation system is 0.67 mol/L; a spherical micelle is thus constructed. The paper is organized as follows: A brief description of the model potentials is given. Details of these potentials have been discussed in ref 1. The construction of the model system and the simulation run are outlined. The time evolution of the shape and of the degree of neutralization of the micelle surface charge by counterions is discussed. The structure is analyzed in terms of density and probability functions with respect to the center of mass of the micelle. From these functions an interface thickness, a surface roughness, and a surface area per head group on the surface are determined. A nearest-neighbor density function for the chloride ion-nitrogen atom pair is calculated to describe the association of counterions with the micelle surface. The conformation of the alkyl chains is analyzed in terms of trans-gauche ratios. The orientation of alkyl chain bonds is determined in terms of a bond order parameter. 0 1994 American Chemical Society
A n-Decyltrimethylammonium Chloride Micelle in Water
Model Potentials The interaction potentials between the particles in the model system are those which have been used in the simulation of a C16TACmonolayer.' The methylene and methyl groups are described as united atoms. The chains are flexible, with intramolecular potentials as functions of torsional angles, bond angles, and bond stretches. The head groups, which have a tetrahedral structure, are rigid. The charge distribution of the head group was obtained by ab initio calculations.' The nitrogen atom carries a charge of -0.56e. The three methyl groups and the methylene group next to the nitrogen carry a charge of 0.39e each. For the intermolecular interactions we employ effective pair potentials. Combining rules are used for the Lennard-Jones terms. The TIP4P water model established by JorgensenZS is used. For thechloridewater interactionswe employ the potentials introduced by Bounds.26 All parameters are explicitly given in ref 1. Model System and Simulation Procedure As mentioned above the model systemconsists of 30 amphiphile molecules, 30 chloride ions, and 2166 water molecules. Theinitial arrangement is based on a configurationof the c16TAcmonolayer simulation:' 15 amphiphiles, 15 chloride ions, and 1083 adjacent water molecules are taken from a selected configuration. This configuration is chosen because of its similarity to a hemispheric part of a micelle. The n-hexadecyl chains are reduced to n-decyl chains. This setup is then duplicated by mirroring at a plane, xy plane, near the air-monolayer interface in such a way that the minimum distance between atoms of the original and the duplicated part is 3 A, roughly the u value of the C-C LennardJones potential. A system of 30 amphiphiles, 30 chloride ions, and 2166 water molecules is thus obtained. This set of particles is surrounded by a periodically replicated simulation cell. The extension of the simulation cell in the z direction is much larger than that in the x and y directions. Therefore, the duplicated part is rotated about the z axis which passes through the center of the xy plane and is moved toward the original part, always keeping the minimum distance of 3 A between atoms of the two halves. A total shift of 4 A is achieved. At this point, the dimensions of the cell are 52 A for the x and y directions and 56 A for the z direction. This procedure keeps the structure of the liquid near to the head groups intact. The initial configuration (0 ps) is shown in Figure la. An external potential
is temporarily applied to the terminal methyl groups in the initial simulation phase in order to obtain a spherical micelle. Therein ro is the center of the simulation cell and k, is a force constant. Test simulationshave shown that no external potentialsare needed for the water molecules. The water phase remains intact during the entire equilibration procedure. The water molecules follow the amphiphiles through their interactions with the head groups. The length of the cell is finally reduced to a value of 42 A in each direction in the first 4 ps of the equilibration run. A value of 5 X J for k, is high enough to move all terminal methyl groups to the interior of the micelle within these 4 ps. The final value for the dimension of the simulation cell is chosen to obtain a normal water density outside of the micelle. The concentration c of CloTAC in the simulated system then amounts to:
The Journal of Physical Chemistry, Vo1. 98, No. 2, 1994 713
interval a time step of 1.5 fs is used, whereas the time step for the other simulation times is 1.O fs. A nearly spherical micelle emerges during this equilibration run. The external potential, eq 1, is then turned off. In the following 2 ps the velocities are scaled, finishing the equilibration run at a given temperature of T = 294 K. This system is then simulated for 215 ps. The last 120 ps, during which 2400 configurations are stored, is used to analyze the structure of the micelle.
Results and Discussion We study first the time evolution of the micelle shape. The principal moments of inertia 11,Iz, and 4 are calculated in order to describe the shape.27 In Figure 2 the ratios I I / I ~I2/I3, , and I l / h , averages over 1 ps, are shown as a function of the equilibration and simulation times. The values for the first 10 ps are a consequence of the reduction of the simulation cell from 5 2 A X 5 2 A X 5 6 A to its finalvalue of 4 2 A X 4 2 A X 42 A. After about 12 ps all three ratios have reached values of about 1.05, indicating a nearly spherical micelle. This shape remains almost unchanged up to a time of about 65 ps, see Figure lb. Such a geometry is to beexpected becauseoftheextemal potential, see eq 1, which is used during the first 57 ps. From 65 ps onward one small and two large moments of inertia appear. A conformational change is observed at about 234 ps. The averages over the last 120ps for the largest, second largest, and smallest moment are 1.31:1.22: 1. From these values it is evident that the micelle is slightly prolate ellipsoid, see Figure IC. The degree of neutralization of the micelle is defined in the same way as in the study of a reversed ionic micelle by Linse.26 A nearest-neighbor density function for the chloride ion-nitrogen atom pair is calculated. The integral of this function at the first minimum is defined as the degree of neutralization. The first minimum is at 6.6 A (see Figure 6). Figure 3 shows the degree of neutralization as a function of the equilibration and simulation time. The values are again averages over 1 ps. From the initial configurationto about 90 ps the degree of neutralizationdecreases, reaching a value of 25%. After this, there are only fluctuations. This seems to indicate that the micelle system is in a state of equilibrium after about 90 ps. The average degree of neutralization for the last 120 ps is 27%. In other works, only about 8 out of 30 chloride ions are in the direct vicinity of a head group. Figure 4 shows the number density profiles pi of the chain members (C2-C10), ofthe headgroups (Cl, N, and three methyl groups),of the water molecules (oxygen atoms), andof thechloride ions with respect to the center of mass of the micelle. The oxygen density profile indicates that the "interior" of the micelle, a sphere of about 8 A in radius, is completely dry. In the region from 8 to 18 A the water density increases smoothly. Between r = 18 and 21 A a region of normal water density po (PO = 0.0334 l/A3) is observed. We recall that this is used as the criterion for the choice of the size of the simulation cell. A normal water density outside of the micelle is necessary to obtain a reasonable density for the chains. Shelley et a1.15 found in an MD study of a sodium dodecylsulfatemicellea water density about 20% above the normal (and a system pressure of about 6 kbar). The chain density was also about 20% higher than obtained experimentally under atmospheric pressure. The hydrocarbon density for the CloTAC micelle here is roughly 10% lower than in their study. The interface thickness Ao, defined according to Carpenter et al.29 as
& = r ( 0 . 9 ~-~' ()0 . 1 ~ ~ )
where Namphis the number of amphiphile molecules, N A is Avogadro'e constant, and Vboxis the volume of the simulation cell. This c value is about 10 times higher than the cmc value. During the next 53 ps the system is equilibrated with k,values reduced stepwise from 1 X to 1 X 10-2' J. During this time
(3) is 6.0 A, for 10.4 A I r I16.4 A. This value corresponds to two to three water layers in the interface region. Because of the ellipsoidal shape of the micelle, this value is an upper limit for the interfacial width. The surface roughness u30,which is defined by a Gaussian approximation to the nitrogen number density
714
The Journal of Physical Chemistry, Vol. 98, No. 2, 1994
B8cker et al.
a
b
C
Figure 1. Instantaneous configurationsof the micelle, omittingthe water moleculesand the chloride ions for visual clarity, as viewed from two mutually perpendicular directions at (a) 0 ps (initial configuration), left side, projection onto the xy plane; right side, projection onto the xz plane (b ) 59 ps (end of equilibration), and (c) 275 ps (final configuration).
A n-Decyltrimethylammonium Chloride Micelle in Water .a Q
The Journal of Physical Chemistry, Vol. 98, No. 2, 1994 715
1.8
O
o
5
r; 8 .
0.01
/?
/
0.000 0
50
100
150
200
250
/ ps Figure 2. Ratios of principal moments of inertia of the micelle as a function of the equilibrationand simulationtime: 11/13 (solid line),12/13 (dotted line), and 11/12 (dashed line). t
5
, 4,
;;I
10
,\,
15
22, , ~, 2 0.-.
0
,-A 25
r/A
*
30
Figure 4. Number density profilespi for chain members (C2-ClO) (solid line), for head groups (Cl,N, and the three methyl groups) (dotted line), for water molecules (oxygen atoms) (dashed line), and for chloride ions (chain dotted line) with respect to the center of mass of the micelle.
c
0
Figure 3. Degree of neutralization of the micelle as a function of the equilibration and simulation time.
function p ( r ) ,
has a value of 2.5 A. The head groups of the micelle are thus located in a molecularly sharp interfacial region. This result is in very good agreement with a neutron study by Bradley et al.31 of an aqueous solution of CloTAB. These authors found that a smooth layer with a thickness between 2 and 3 A is formed if the area per head group is larger than about 65 A2 on a monolayer or on a micelle surface. Other authors,3" investigating n-alkyltrimethylammonium micelles with different chain lengths in aqueous sodium bromide or sodium chloride solution, also found values between 2 and 3 A. Comparing the oxygen density with the head group density profiles shows that the water molecules penetrate into the head group region, but not into the region of the chains. The radial overlap between oxygen atoms and chain members is caused by the location of some chain members in the head group region. From the number density functions pi probability distributions Pi are calculated with the following equation:
Pi(r) dr = 4?r2pi(r) d r (5) Figure 5 shows the probability distributions of C10, C5, N , and C1-. A mean distance ( m )from the micellar center of mass of 14.2 A is found for the nitrogen atoms. Assuming a spherical surface this leads to an area per head group of 84 142. According to a study of micelles by Imaeet a1.8 the aggregation number Namounts to 44 ( C I ~ T A C )62 , ( C I ~ T A C and ) , 84 (c16TAC). The area A per head group is reported as 85 A2 (C12-
5
10
15
20
25
D
r/A Figure 5. Probability distributions Pi for carbon atoms C10 (solid line), for carbon atoms C5 (dotted line), for nitrogen atoms N (dashed line), and for chloride ions C1-(chain dotted line) with respect to the center of mass of the micelle.
TAC), 80 A2 (CldTAC), and 78 A2 (ClbTAC), and the mean radius (r) of a micelle is 17.1 A (ClZTAC), 19.9 A (CI~TAC), and 22.6 A (C16TAC). The mean radiusof a micelle is determined with an error af about hO.5 A. These values are estimated under the assumption of a spherical micelle. Berr et alO7observed an increase of the aggregation number with increasing concentration of amphiphile molecules. Comparing these two studies738 and a study of Brady et al.32 we conclude that these aggregation numbers can be determined with an uncertainty of about f 5 amphiphile molecules. The assumed aggregation number of 30, the area per head group of 84 A2,and the mean radius ( r ~ of) 14.2 A, observed in our study, are thus in good agreement with the study by Imae et a1.,8 taking into account the different chain lengths. The widths of the distributions Pi of the nitrogen atoms and of the Cl-CS methylene groups are nearly identical, whereas the distribution of the C10 groups is much broader, indicating that the flexibility strongly increases toward the end of a chain. The curves in Figure 5 agree essentially with the traditional picture of a micelle, where the head groups are distributed near the surface of a sphere and the hydrophobic core shows liquid-like character. In Table 1 the average radial position (ri)is given for all chain members, the nitrogen atoms, and the chloride ions. The difference of the mean radii between nitrogen atoms and chloride ions is 6.9 A. This fact and the comparison of the probability distributions (Figure 5 ) suggest that if a chloride ion is associated to a head group it is located outside of the micelle. Figure 6 shows the nearest-neighbor density function for the chloride ion-nitrogen atom pair and the integral of this function. Note that a value of 1 in the integral corresponds to 100%
B&ker et al.
716 The Journal of Physical Chemistry, Vol. 98, No. 2, 1994
P
10 0
80-
-06
60-
2 \
9
40-
00
,
o
Figure 8. Orientation of a bond relative to a radius vector r. 0 is the center of mass of the micelle.
-02
20-
,
2
,
, ,
, ,
4
,
1
e
,
,
I
.
,
,
I
.
,
,
I
,
'
,
iz
io
e
,
14
'
, ,
18
Figure 6. Nearest-neighbor density function CI--N. The jagged line gives the density function, and the smooth line gives the integral of the density function.
\
\ ......
.............................................................
-0.1
11
2
3
4
e bond number 5
7
8
s
J
10
Figure 9. Bond order parameter S as a function of individual bonds in a chain.
0
1
2
3
4
5
6
7
8
number trans conformations Figure 7. Probability of finding a given number of trans bonds in a chain.
TABLE 1: Average Radial Position (q) for CI; N, and Each Chain Member, with Respect to the Center of Mass of the Micelle ~~~~~~
i
ClN
c1
c2 c3 c4
(ri)/A 21.1 14.2 13.2 12.2 11.3 10.6
i
c5 C6 c7 C8 c9 c10
(ri)/A 9.9 9.2 8.6 8.2 7.8 7.6
neutralization. The bimodal form of the density function is quite obvious. The first minimum is at r = 6.6A. The integral at this radius gives a degree of neutralization of 27% as discussed above. It is pointed out that the correspondingvalue for the monolayer' is 74%. This remarkable difference is clearly a result of the larger area per head group on the micelle surface in comparison to the monolayer. The chloride ions in the micelle system are less attracted because of the low charge density. The weak second maximum in the density function implies that there is virtually no second solvation shell associated with the micelle. Assuming a radius of 9.3 A, which was found in the monolayer' analysis, 37% of the chloride ions are found in the region between 6.6 and 9.3A, leaving 36% free chloride ions. In the first solvation shell, 93% of the chloride ions have one and 7% have twonitrogen neighbors within their shell. Thus, a chloride ion in the first solvation shell is almost always associated with only one amphiphile molecule, in contrastto the monolayer' result. The trans-gauche ratio is calculated in order to describe the conformation of the chains. A bond is considered to be in a trans stateifcosg
