Langmuir 1992,8, 3155-3160
3155
Molecular Dynamics 1nvesti.gationof a Newton Black Film Zulema Gamba,+Joseph Hautman, John C. Shelley, and Michael L. Klein* Department of Chemistry, The University of Pennsylvania, Philadelphia, Pennsylvania 19104-6323 Received June 1, 1992. In Final Form: September 28, 1992
Soap films consisting of amphiphilicmolecules and water have been investigatedby molecular dynamics simulations. The model systems studied contained a specified number of water molecules sandwiched between two monolayers of sodium dodecyl sulfate (SDS) molecules whose alkyl chains were confined to a fixed in-plane densitythrough the imposition of periodicboundary conditions. The constituentsinteracted via potentials taken from a previous study of an aqueous SDS micelle. The calculated overall width and the variation of the density profile across the film are consistent with deductions based on recent X-ray reflectivity data. Equilibration from an ordered initial (far from equilibrium)structure involves extensive conformational changes in the alkyl chains, molecular displacements, and rearrangement of the water molecules around and behind the -0-SOa headgroups. Films of different thickness (i.e. differing degrees of hydration) were stable over the time scale of our simulationswhich spanned hundreds of picoseconds. On the basis of energetic considerations,and comparison to experiment, the Newton black film studied in the X-ray work likely contained from four to six water molecules per surfactant. Introduction
Bilayers of soap molecules, formed by drawing a film from an amphiphile/water solution, thin to stable free standing structures with aspect ratios, diameter to film thickness, on the order of lo6 or more. In the final stage of thinning the films do not reflect light and are referred to as “blackfilms”.l Recent X-ray reflectivity experiments measured the thickness of black films formed from sodium dodecyl sulfate (SDS) molecules and water.2 Two types of black films were observed. The first, identified as a “common” black film (CBF), was deduced to be about 54 A thick. However, in CBF films the equilibrium thickness is known to exhibit a strong dependence on the salt concentration in the soap solution and ranges from the value reported in ref 2 to -2000 A in very dilute systems.l The second film, a “Newton” black film, was measured to be about 33 A thick,2and this type of film appears to have a minimal dependence on salt c0ncentration.l The Newton black film (NBF) is generally found only in the presence of higher salt concentrations and is thought to represent the maximally thinned bilayer structure. In either the CBF or NBF case, the measured thickness suggests that there is room for only a few layers of water molecules separating the two opposing monolayer^.^,^ A schematic diagram of the cross section of a NBF is shown in Figure 1.
We have recently developed a potential model to describe SDS [Na+CH3(CH2)110S03-] molecules in aqueous s~lution.~ This “forcefield”was employed in molecular dynamics simulations on both a monomer and a micellar aggregate with some success. In the present article, we have used the same potentials to investigate Newton black films, again, with molecular dynamics simulations. Anticipating our results, we will see that the model yields black film bilayers of SDS which are stable on the 100 ps time scale of the simulation. Comparison of the measured density profile with simulation results for systems with + Permanent address: Departamento de Fisica,ComisionNacional de Energia Athmica, Av. Libertador 8250, (1429) Buenos Aires, Argentina. (1) Jones, M. N.; Mysels, K. J.; Scholten, P. C. J. Chem. SOC.,Faraday Trans. 1966,62, 1336. ( 2 ) BBlorgey, 0.;Benattar, J. J. Phys. Reu. Lett. 1991, 66, 313. (3) DeFeijter, J. A.; Vrij, A. J. Colloid Interface Sci. 1979, 70, 456. (4) Shelley, J.; Watanabe, K.; Klein, M. L. Int. J. Quantum Chem: Quantum Biol. Symp. 1990,17,103.
Figure 1. Schematic diagram showing, in profile, the gross structure of a black film. The diagram contains distinct regions consisting of (a) water molecules and counterions, (b) chain headgroups, and (c) hydrocarbon tails.
varying degrees of hydration suggests that the Newton black film studied in the X-ray scattering experiment contained between four and six water molecules per surfactant. The structure of the film is examined in detail with particular emphasis on the conformations of the alkyl chains and the nature of the water in the core of the black film. The headgroup hydration and the counterion condensation are also discussed. Method Molecular dynamics calculations for an idealized Newton black film were performed in the (N,V,T) ensemble at T = 298 K.5 The systems studied consisted of either 36 or 64 dodecyl sulfate anions, an equal number of Na+ cations, and a number of water molecules, n,, ranging from zero to six per SDS molecule. As mentioned above, the Newton black film is stable in the laboratory only when formed from a solution containing a sufficient concentration of a salt such as NaC1.l Nevertheless, the thinned film is thought to contain only the SDS molecules and water but no salt.3 The starting configuration for the calculations was constructed by arranging half of the chain molecules on a body centered square lattice of 3 X 3 or 4 X 4 unit cells, in the (x,y) plane, with lattice constant a = 8.124 A. Molecular centers of mass of the SDS anions were initially placed at (0, 0, 2,) and (l/2, 1/2, zc), respectively, which resulted in a surface density, 33 A2 per alkyl chain, corresponding to experimental estimates.2 The head group (5) Allen, M. P.; Tildesley, D. Clarendon: Oxford, 1987.
J. Computer Simulation of
0743-7463/92/2408-3155$03.00/0 0 1992 American Chemical Society
Liquids;
3156 Langmuir, Vol. 8, No. 12, 1992
lattice in the starting configuration defined a surface in the (x,y) plane and the initially all-trans hydrocarbon chains were tilted about 55O from the positive z direction. A second layer, which completed the black-film bilayer, was obtained by reflection in the z = 0 plane. Headgroups f7 A, which resulted were located in two planes at z in an initial distance between the SDS layers about twice that deduced from the X-ray measurements. The Na+ counterions were placed in the plane z = 0, using the same square lattice but with a basis (1/4, 1/4, 0), (1/4, 3/4, 0), (3/4, l/4,0) and (3/4, 3/4, 0). The water molecules were initially distributed in a regular array between the headgroups. This starting configuration is far from equilibrium, but in only a few picoseconds of calculation the system compresses spontaneously to a film that, after thermalization, is stable for simulation runs of several hundreds of picoseconds. The equationsof motion were integrated using the Verlet s . ~Periodic algorithm with a time step of 3.5 X boundary conditions were used in all directions. The side of the simulation cell in the direction perpendicular to the bilayer plane ( z ) was 100 A, i.e. about 3 times the dimensions in the ( x , y ) plane of the bilayer and about 3 times the final film thickness. The electrostatic interactions were handled via the Ewald method, and sums in real space are performed with a 12-A~utoff.~ The LennardJones atom-atom parameters were the same as those used in a simulation of a quasi-spherical SDS micelle in aqueous s ~ l u t i o n .The ~ chain molecules consisted of negatively charged -&so3 headgroups attached to a single chain of 11methylene groups terminated by a methyl group. The methylene and methyl groups were represented by spherical pseudoatoms with rigid C-C bonds constrained using the Shake a l g ~ r i t h m . The ~ , ~ molecule was otherwise fully flexible with bond bending and torsion potentials. Reference 4 contains full details of the intramolecular and intermolecular potential parameters for the water, Naf counterions, and surfactant molecules. Before choosing the number of water molecules per chain, n,, to be used in a large simulation with 64 SDS molecules, we performed a series of calculations using a smaller system of 36 surfactants, and n, ranging from 0 to 6. From these calculations we estimated that the system should contain around six water molecules per SDS in order to reproduce the experimental thickness of the film. Typically, the equilibration period was about 50 ps, with most of the rearrangement of molecules occurring during the first 5 or 10 ps. While the overall chain tilt developed rapidly, the conformational equilibrium of the alkyl chains is attained more slowly. The evolution of the gauche defect concentration for representative torsional angles is shown in Figure 2. Atomic coordinates were usually stored for subsequent analysis. The total length of the run with the larger system exceeded 200 ps.
-
Results A. Film Profile. All of the systems studied preserved the gross bilayer structure, similar to that shown schematically in Figure 1, for the entire length of the stimulation. A snapshot from the simulation with 64 chains and n, = 6 is shown in Figure 3. Although the average structure exhibits well-defined water, headgroup, and alkane chain regions, there is considerable roughening in the various interfaces. In particular, Figure 1 is too idealized since water molecules envelop the headgroups, forming a well-defined solvent sheath. The Na+ coun-
Gamba et al.
,, ,
1
I
( (
,I
,,
" ( 1
ly""
"'1
7 4
C
2 .8 m L
4
.6
1
t
.4' 0
1 '
'
'
' 20
'
'
'
I
40
'
'
'
I
60
'
'
t(ps.)
'
' 80
Figure 2. Evolution of selected alkyl chain torsional angle (7) populations over a 80-ps segment of a molecular dynamics trajectory for a sodium dodecyl sulfate black film (see text). The dihedral angles are numbered starting from the headgroup; r4is the first unconstrained angle.
terions are found both as solvated ions in the water core region and as associated counterions (contact-ion pairs) in the headgroup region. These structural features are similar to local configurations observed in the SDS micelle ~imulation.~ Although there are many gauche defects present in the alkyl chains, a clearly discernible average tilt is present (about 45O) which is consistent with the chain tilt versus area-per-chain data for long-chain n-alkane monolayers at somewhat higher in-plane densities.' Analysis of X-ray reflectivity measurements of ref 2 yields the electron density profile which is compared to the calculated profile in Figure 4. As suggested by the small system studies, the overall width of the bilayer is well reproduced with n, = 6. The experimental analysis separated the densityprofile into three regions, as in Figure 1, where (a) corresponded to water/Na+ (b) corresponded to headgroups, and (c) corresponded to hydrocarbon tails. The molecular dynamics results confirm that this basic structure is reasonable but the calculated widths of the separate regions differ somewhat from experiment, however, this point will be discussed later. Figure 5 shows selected atomic density profiles of the water molecules, Na+counterions, headgroups, and tails. The water density profile clearly overlaps the chain density-another indication that the headgroups are hydrated, i.e. at least partially surrounded by solvating water molecules. Although the overall width of the electronic density profile matches the experimental profile, some discrepancies can be noted in Figures 4 and 5. In particular, the width of the water region, bounded by the headgroupdensity peaks, is somewhat wider in the simulation while the hydrocarbon region is narrower. The discrepancy in the core region of the film suggests that there should be less water between the headgroups or that the calculated roughness of the layer formed by the headgroups is much larger than in experiment. In fact, as can be seen in Figure 6, this central (core) region of the black film is better fitted by the profile from the simulation with n, = 4. The fact that the hydrocarbon region is too narrow is consistent with the calculated average chain tilt which, in turn, follows from the assumed surface density of alkyl chains of the SDS molecules. A slightly thicker hydrocarbon region would likely result from an all-atoms description of the alkyl hai ins.^^^ B. Chain Mobility and Conformational Equilibrium. The chain molecules in the NBFs are known to be
(6) Ryckaert, J. P.; Ciccotti, G.;Berendsen, H. J. C. J.Comput. Phys. 1977, 23, 321.
(7) Bareman, J. P.; Klein, M. L. J. Phys. Chem. 1990, 94, 5202.
Langmuir, Vol. 8, No. 12,1992 3157
Molecular Dynamics Investigation of a Newton Black Film
t
F
Figure 3. Orthogonal views of a typical configuration of the black film taken from the molecular dynamics calculation with 64 SDS surfactants and 384 water molecules (not shown). The Na+ ions (green), 0 atoms (red), S atoms (yellow), CH2 groups (blue), and terminal methyl groups (purple) are drawn with covalent radii.
mobile? However,the diffusion coefficient is rather small; the measured value for SDS is on the order of 0.06 X
cm2s-1. Diffusion at this rate would be difficult to observe in the present calculations because, on average, each
3158 Langmuir, Vol. 8, No. 12,1992
Gamba et al.
.8
.6
N v
.4
n .2 I
I 1 ' I '
0 0
I "
I
10
"
"
I
20
'
"
" " '~ 30, ( A )
40
Figure 4. Electron density profile (electrons/A3)derived from the simulation of a black film having 64 SDS molecules and 6 water molecules per surfactant. The dashed line is the density profile derived from the X-ray reflectivity data of ref 2. The overall width of the film is in good accord with experiment but the central core region is too wide and the hydrocarbon region too narrow (see text).
20
10
0
30,
( A )
40
Figure 6. Electron density profile (electrons/A3)derived from a simulation of a black film having 36 SDS molecules and 4 water molecules per surfactant. The dashed line is the density profile derived from the X-ray reflectivity data of ref 2. Overall this film is too narrow, primarily because of the width of the alkyl chain region (c) (see text).
.04
,.03 N v
.02
.o 1 4
6
10
8
12
14
Tn 0 0
10
20
3 0 2 ( A ) 40
Figure 5. Atomic density profiles for water molecule 0 atoms (dashed line), Na+ counterions (dotted line), S atoms of the headgroups (full line), and terminal methyl groups (long-dashed line), derived from a simulation of a black film having 64 SDS molecules and 6 water molecules per surfactant. The degree of water penetration into the headgroups region should be noted. The vertical scale refers to the water density in units of A-3.Note that the value at the center is close to that of bulk water (0.033 A-3). The values of the other density profiles have been scaled by a factor of 6 relative to the water profile.
headgroup should diffuse only about 5 A in a 1-ns calculation! This leaves open the question as to the degree of equilibration of the black film systems studied herein whose molecular dynamics calculations were run for times of only a few hundred picoseconds. Although the requisite nanosecond run times are feasible at the present time, the expense will likely not be justified until the experimental data and the interaction potentials have been further refined. Another potential difficulty in the simulation of the black film is the expected long time scales required to achieveconformational equilibration of the SDS alkyl tails, a problem common to all chain systems.lO We have already alluded to the fact that frequent transitions occur between conformers (recall Figure 2). In the present case we have (8) Bareman, J. P.; Klein, M. L. In Interface Dynamics and Growth; Materials Research Society Symposium Proceedings; Materials Research Society: New York, 1991; Vol. 237 p 271. (9) Clark, D. C.;Dann, R.; Mackie, A. R.; Mingins, J.; Pinder, A. C.; Purdy,P. W.; Russel1,E. J.;Smith,L. J.;Wilson, D. R.J. Colloid Interface Sci. 1990, 138, 195. (10) de Loot, H.; Harvey, S. C.; Segrest, J. P.; Pastore, R. W. Biochemtstry 1989, 30, 2099.
Figure 7. Fraction of gauche g+(circles) and g7triangles) conformers as a function of, T", the dihedral angle number. The first three angles, involvingthe headgroup atoms are constrained by steric effects to be either trans or gauche. Values are shown for two separate 42-ps segments, labeled (a) and (b), taken from the molecular dynamics trajectory for the system shown in Figure 3.
monitored the fraction of gauche g+ and g- conformers throughout the simulation. These values are shown in Figure 7 for two 42-ps segments taken from the run on the 64 surfactant system. Fhe g+ and g- fractions appear to have stabilized fairly well. It can be seen in Figure 3 that although the headgroups are disordered, in orientation and location, the same is not true of the chains which have roughly the same overall orientation and tilt angle. This fact is consistent with studies on other chain systems at roughly the same density.' C. Headgroup Hydration. Figure 8 includes several pair correlation functions g(Atom-Atom) relevant to the study of the local structure of the film. The function g(SS) shown as the dashed curve has a broad first peak centered at - 6 A, while that of g(Na+-Na+), shown as a dotted curve,has a more distinct peak at -4 A. In contrast, the g(S-Na+) function, shown as a full curve, has a very sharp bimodal first peak at small distances. Since there are 1.7 neighbors in the whole (double) peak, surfactant headgroups must share counterions. Analysis of the orientational probability distribution of Na+ around the headgroup indicates that the -3-A peak in Figure 8 is due to counterions behind the headgroup and that the next peak is mainly due to counterions on the terminal 0 atoms. The same analysis can be carried out for the water around the headgroups. The pair distribution functio-; for S-Ow,
-
Langmuir, Vol. 8, No. 12, 1992 3159
Molecular Dynamics Investigation of a Newton Black Film 2.5 -7.7
2
h
p.
-
1.5
v
J
f
-8
1
-8.1
"
0
'
I
"
'
2
"
4
' 6
n,
2
0
8 r
6
4
(A)
10
Figure 8. Pair distribution functions g(Atom-Atom) for S-S (dashedline), Na+-Na+ (dotted line) and S-Na+ (full line) atom pairs for the system in Figure 3. 32
u
24 0 22
2
4
6
n,
Figure 9. Overall f i i thickness,taken tobe the distance between the peaks in the terminal methyl group distributions (see Figure 5) versw, h, the number of water molecules per surfactant based on the present molecular dynamics calculations. The dashed line shows the expected asymptotic result correspondingto the bulk water model employed in the present simulations (ref 13). The solid circles are from n, = 36 and the open circle from n, = 64.
5-H,, and their angular distributions confirm that there are indeed some water molecules behind the heads (recall Figure 3). D. Water Mobility. The calculated mobility of the water in the Newton black film is much less than that of bulk water at room temperature. Typical calculated values of the average diffusion coefficient for the water molecules are an order of magnitude lower than for bulk water. It should be noted, however, that the environment for the water molecules is quite inhomogeneous. A calculation of the diffusion coefficient, as a function of the distance to the central plane, shows that the mobility of the water is about half the bulk value in the middle of the black film. This finding is consistent with previous calculations on concentrated soap solutions.ll For thicker water layers coated by surfactant (Langmuir) films, the calculated diffusion coefficient in the middle of the central water region is comparable to that of bulk E. Equilibrium Water Content. As a fiist step toward determining the number of water molecules preferred in the model system, we have undertaken a study of the energetics of the bilayers as a function of the number of water molecules per chain, n,. This investigation of film structure and energy was carried out using the smaller system with N c k = 36. Surprisingly, a comparison of results from runs with Nchain= 36 and 64, for n, = 6, appears to show only minor effects due to the system size. (11) Watanabe, K.; Klein, M. L. J . Phys. Chem. 1991, 95, 4158. (12)Raghevan, K.; Rami Reddy, M.; Berkowitz, M. L. Langmuir 1992, 8, 233.
Figure LO. Plot of potential energy of the black film per SDS molecule, u, versus h,the number of water molecules per surfactant in the black film. The dashed curve is a quadratic best fit to the data.
In particular, the structure of the hydrocarbon region of the black film appeared to be nearly independent of b, and variation in the number of water molecules resulted in a nearly linear increase in overall film thickness (Figure 9). Thus the structure of the hydrocarbon region of the black film is largely independent of the water content. Knowledge of free energy differences for the varioue systems is required in order to determine the preferred hydration state of the black film. Experiments show that these films are indeed in thermodynamic equilibrium.' This implies that the water between the layers is in equilibrium with a reservoir of bulk water at the edge of the bilayer region. The free energy for the composite bilayer-reservoir system containing a total of N water molecules is then a sum of the bilayer free energy CANw) and the reservoir free energy (N-Nw)h,where h is the chemical potential of the bulk water and N , is the number of water molecules in the bilayer
G,t(N,) = GfW,) + (N-N,)PO (1) Subtracting the constant N h , we obtain an energy In the thick film limit of large N,, G(N,) corresponds to twice the surface excess free energy for a monolayer film
-
GW,) 2G, (3) where G, is the free energy difference between the monolayer water system and the "bulk" energyof the water. Simulations in the grand canonical ensemble would provide a natural tool for understanding the subtle energetics underlying the formation of NBFs and CBFs. For the present, first study we have not calculated free energies but only the potential energy differences for the bilayer system as a function of n, = N w / N c k .However, with crude assumptions about the n, dependence of entropic contributions to the free energy, G, this is sufficient to gain insight into the factors determining the equilibrium thickness of the black film. In Figure 10 we have plotted the potential energy component u(n,), of the free energy defined by eq 2 which is normalized by the number of SDS molecules u(nJ = u@,) - nwuO (4) where n, is the number of water molecules per SDS and uo is the potential energy of a water molecule in the bulk (UO = -46.6 kJ/mol for the SPC/E model's). If entropic contributions to the free energy could be neglected, this (13) Berendsen, H. J. C.; Grigera, J. R.; Staatama, T. P. J. Phys. Chem. 1987,91,6269.
Gamba et ai.
3160 Langmuir, Vol. 8,No. 12,1992 trend would predict a local equilibrium value of n,around 3.6 1. Since we observe the structure of the chains in the bilayer to be largely independent of water content, the relevant entropy contribution is then the difference between the entropy of water in the bilayer and that in the bulk. For thick black films this difference clearly approaches zero. For the thinnest films the water in the core is frozen (or glassy) so that the entropy difference will tend to make the free energy for smaller values of n, larger. The equilibrium hydration value would thus be shifted to a value of n,larger than that predicted on the basis of u alone. In all of the present simulations,however, the mobility of the water is far less than that of bulk water and it is likely that the entropic contribution does not vary appreciably.
*
Discussion The present model has been used to study the structure of a sodium dodecyl sulfate Newton black film for times short compared with the characteristic diffusion time for surfactant monomers. The molecular dynamicssimulation qualitatively reproduces the canonical bilayer structure of the black film with distinct hydrocarbon tail, headgroup, and water regions.3 The simulation indicates, however, a significant roughness to the interfaces between these regions reflecting both disorder in the headgroup layer and penetration of water molecules into the headgroup region. Quantitative discrepancies between the simulation results and experimental data have already been noted. These may be due to inadequacies in the potential model: exceedingly long conformationalequilibration times, and/ or uncertainties in the value of the surface area per chain.2 In the experimental analysis of ref 2 the area per chain is deduced from the width of the hydrocarbon region and the measured electron density, Pe = 0.27e/&. The latter value is significantly lower than densities reported in Langmuir monolayers of lon hain carboxylic acid mole c d e ~namely, , Pe 15 0.33 e/ &4 and the value obtained in the present simulation, 15 0.32e/A3. The experiments of ref 14 and studies using similar potential modeling in long-chain monolayers7show that changing the area per (14) Kjaer, K.; Ale-Nielwn, J.; Helm,C. A.; Tippman-Krayer, P.; Mbhwald, H.J. Phys. Chem. 1989,93,3200.
chain does not, to fiist order, affect the density in the hydrocarbon region but simply changes the tilt angle of the chain tails. We might therefore expect that any modification of the surface area per surfactant in the present simulation would lead to different widths of the hydrocarbon region (b) in the density profilee7 The reported low value for the hydrocarbon density could perhaps be obtained in the present model if a much higher number of conformational defecta were to form over a time scale much longer than a nanosecond. Actually, the experimental density profile can be integrated to obtain the total number of electrons in the black film. If one subtracta the contribution from the chain molecules and the Na+ counterions and attributes the remainder to water, one obtains n, = 2.6 molecules per SDS, which is consistent with the simplistic estimates based on the analysis of the potential energy given above. The various estimates of the degree of hydration present in the SDS Newton black film bilayer given herein seem to be a variance with the commonly held view16 that the minimum level of hydration of SDS is rt, = 10. The situation in micellar solutions may perhaps be somewhat different from that in the Newton black film. Finally, the results obtained in the present simulation studies of a black film bilayer are sufficientlyencouraging that it would likely be worthwhile to attempt a more detailed and thorough study in order to probe the nature of the hydration forces between the surfactant bilayers.16 Such a study would likely benefit from the use of an efficient algorithm for sampling the conformations of the alkyl chains.17
Acknowledgment. We are grateful to J.-J. Benattar, Bob Thomas, and Helmuth M6hwald for their comments and suggestions. The research described herein was supported by the National Institutes of Health under GM 40712 and RR04882. Some of the calculationswere carried out at CNSF. Registry No. SDS,151-21-3. (15) Bleasdale, T. A.;Tiddy,G. J. T. In The StructureandEquik'bn'um Properties of Colloidal Systems; Kluwer: Dordrecht, 1990; pp 397-414. (16) Ieraelachvili, J. N.; Wennerstrt", H.Langmuir 1990,6, 873. (17) Siepmann, J. I,; Frenkel, D. Mol. phy8. 1992, 76, 69.
