Langmuir 2004, 20, 5127-5137
5127
Computer Simulation Studies of Newton Black Films Fernando Bresme* and Jordi Faraudo† Department of Chemistry, Imperial College London, Exhibition Road, London, SW7 2AZ, United Kingdom Received October 28, 2003. In Final Form: March 25, 2004 We study via molecular dynamics simulations thin films (Newton black films, NBF) consisting of water coated with sodium dodecyl sulfate (SDS) surfactants. We analyze in detail the film properties (distribution of particles, pair correlation functions, roughness of the film, tilt angle of the hydrocarbon chain, electron density profiles, and mobility of water molecules) as a function of water content in the film core (i.e., film thickness, H). Our simulations indicate that water is part of the bilayer structure as solvation water. We estimate that around 2.25 water molecules per sufactant are part of this solvation structure. The structural analysis of the NBF shows that the headgroups exhibit a high degree of in-plane ordering. We find evidence for the existence of cavities in the monolayer, where only water is present. The basic structure of the monolayer is conserved down to water contents of the order of 4 water molecules per surfactant (H ≈ 11 Å). The computed monolayer roughness for the present model is 2.5 Å, in good agreement with the experimental data. We find that the roughness is very sensitive to the details of the interatomic potentials. Water mobility calculations emphasize the sluggish dynamics of very thin NBF. Diffusion coefficients of water in the lateral direction strongly decrease with film thickness. We find that the typical mean squared displacement of water in the direction normal to the bilayer is between 9 and 80 Å2. Overall, our results indicate that the equilibrium SDS Newton black films studied in the X-ray experiments contain from 2 to 4 water molecules per surfactant.
I. Introduction Black films stabilized by ionic surfactant have been the subject of scrutiny in many experimental and theoretical studies.1-8 These surfactant bilayer structures constitute the basic building blocks of foams, emulsions, and also biological membranes. Black films play an important role as experimental systems where theories of intersurface forces can be tested.9,10 These systems are particularly suited to study questions of relevance in biochemistry and biophysics, such as membrane stability. It is important to point out that the main physics appearing in biological membranes is also present in black films; this makes them useful tools for one to understand the behavior of freestanding membranes, for instance. Also, black films are important from an industrial perspective, because their stability determines the stability of foams and microemulsions.11 The thinnest equilibrium states of soap films were first observed by Hooke and Newton. Depending on their thickness, they are called common or first (CBF) and Newton or second black films (NBF). Typically, CBFs have a thickness of the order of 10 nm, whereas NBFs are * To whom correspondence should be addressed. E-mail:
[email protected]. † E-mail:
[email protected]. (1) Mysels, K. J. J. Phys. Chem. 1964, 68, 3441. (2) Clunie, J. S.; Corkill, J. M.; Goodman, J. F. Faraday Discuss. 1966, 42, 34. (3) Jones, M. N.; Mysels, K. J.; Scholten, P. C. Trans. Faraday Soc. 1966, 62, 1336. (4) Be´lorgey, O.; Benatar, J. J. Phys. Rev. Lett. 1991, 66, 313. (5) Benatar, J. J.; Daillant, J.; Berlorgey, O.; Bosio, L. Physica A 1991, 172, 225. (6) Sentenac, D.; Schalchli, A.; Nedyalkov, N.; Benatar, J. J. Faraday Discuss. 1996, 104, 345. (7) Sentenac, D.; Dean, D. S. J. Colloid Interface Sci. 1997, 196, 35. (8) Bergeron, V. J. Phys.: Condens. Matter 1999, 11, 215. (9) Israelachvili, J. Intermolecular and Surface Forces, 2nd ed.; Academic Press: New York, 1991. (10) Evans, D.; Wennerstro¨m, H. The Colloidal Domain, 2nd ed.; Wiley-VCH: New York, 1999. (11) Boek, E. S.; Jusufi, A.; Lo¨wen, H.; Maitland, G. C. J. Phys.: Condens. Matter 2002, 14, 9413.
roughly 4 nm thick.4 Recent experiments12 showed that the CBF to NBF transition is a first-order phase transition. A challenging problem concerning these systems is the prediction of the equilibrium thickness and the identification of the forces determining such thickness. In the case of CBF, the equilibrium thickness is interpreted as a balance between van der Waals attraction and electrostatic double-layer repulsion forces, as expected from DLVO theory.8 In the case of Newton black films, their stability is not well understood. For very small thicknesses, where the distance between the two monolayers in the black film is very small, new forces not accounted for by the DLVO theory do appear. It is believed that repulsive entropic forces (due to the confinement of protusions of the adsorbed layers) might be responsible for the stability of the film,8,10,13 counteracting the attractive van der Waals forces. Also, the DLVO description of electrostatic interactions (based on the Poisson-Boltzmann equation) has serious limitations at small distances between the charged interfaces.10,13 At small distances, the discrete nature of the monolayer charges becomes important. This effect can be included in the classical DLVO theory.14 Theoretical calculations beyond the DLVO theory and computer simulations7,10,15 show that the strength of electrostatic interactions can be less repulsive than the predictions of the Poisson-Boltzmann equation. However, these results apply to the case of electrostatic interactions between charged rigid surfaces. It is important to remark that interactions between soft interfaces can be very different due to the effect on the charge distribution (and consequently the electric field) of protusions, undulations, and other phenomena typical of soft interfaces. To the best of our knowledge, no clear and satisfactory answer has been provided to the question of how and why NBF are stable. The forces operating in bilayer thin films can only be understood if the film structure is well characterized. For (12) Casteletto, V.; et al. Phys. Rev. Lett. 2003, 90, 048302. (13) Israelachvili, J. Nature 1996, 379, 219. (14) de Feijter, J. A.; Vrij, A. J. Colloid Interface Sci. 1979, 70, 456. (15) Marcelja, S. Nature 1997, 385, 689.
10.1021/la036026w CCC: $27.50 © 2004 American Chemical Society Published on Web 05/13/2004
5128
Langmuir, Vol. 20, No. 12, 2004
this reason, a great deal of experimental effort has been made to understand the structure of Newton black films. The X-ray experiments by Benatar and co-workers4-6 on water films coated by sodium dodecyl sulfate (SDS) have shed light on the complex structure of these systems. The interpretation of the experimental data leads to a NBF structure consisting of a well-defined surfactant bilayer coating a thin water slab approximately 4 Å thick, the total thickness of the NBF being 33 Å. The experiments indicate that the NBF thickness does not change with ionic strength.4 This gives support to a structure in which the only ions present in the film are the counterions, Na+, that preserve the electroneutrality. Interestingly, the experiments also indicate that the aliphatic structure of common black films, with a water slab that is 27 Å thick, and the aliphatic structure of the NBF are essentially the same. The X-ray experiments and diffraction experiments, in general, rely for their interpretation on models that are built a priori using chemical intuition. The validity of these models is tested a posteriori, by comparing the model predictions with the experimental data. Therefore, an independent knowledge of the structure of these complex interfacial systems is desirable. Insight into this structure is needed to understand the nature of the forces operating in thin films at different thicknesses. A route to obtain such structure is computer simulations, which provide a detailed description of the positions of the atoms in terms of well-defined interatomic forces. The first attempt to simulate a Newton black film with the same composition as that considered in Be´lorgey and Benatar’s experiments4 was the work by Gamba et al.16 These simulations predicted interesting structural aspects of the film and pointed out the sluggish dynamics exhibited by NBF. This represents a problem for the molecular dynamics simulations of these systems, and very long runs are needed to obtain equilibrium bilayers. In this paper, we extend further the computer simulation of these interesting systems, taking advantage of modern supercomputing facilities. Specifically, we consider very thin films of water coated with sodium dodecyl sulfate (SDS). We have chosen SDS because this system has been analyzed in several experiments,4,5 and also because SDS is one of the most widely used detergents, with many scientific and industrial applications. The results obtained from this study should nonetheless be of interest to other ionic surfactants, because the physics of the problem should be essentially the same. We find that simulations spanning the nanosecond time scale are necessary to obtain accurate statistics of the equilibrium properties of these complex systems. Our main goal in this paper is to analyze the film properties: pair correlation functions, particle density profiles, tilt angles, electron density profiles, and water mobility for different film thicknesses (different water content). II. Potential Models and Simulation Methods We have performed molecular dynamics simulations of water films coated with SDS surfactants [Na+CH3(CH2)11OSO3 ] in the canonical ensemble at T ) 298 K. As explained in the Introduction, no salt is observed in the NBF film; therefore, the only counterion present is Na+ (salt-free conditions). In our simulations, we explicitly model the water molecules, the Na+ counterions, and the DS- anions (i.e., [CH3(CH2)11 OSO3 ] anions). There are several accurate water models available in the litera(16) Gamba, Z.; Hautman, J.; Shelley, J. C.; Klein, M. L. Langmuir 1992, 8, 3155.
Bresme and Faraudo Table 1. Nonbonding Interactions: Parameters for Lennard-Jones Interactions and Partial Charges (qe Is the Absolute Value of the Electronic Charge) group
(kJ/mol)
σ (Å)
q/qe
ref
Na O(water) H(water) CH3 CH2 CH2 bonded to O(ester) O(ester) S O
0.4179 0.6502 0 0.4577 0.4577 0.4577
2.5840 3.1660 0 3.816 3.816 3.816
1.0 -0.8476 0.4238 0 0 0.137
20 19 19 22 22 22
0.7133 1.0460 0.8786
3.3224 4.0000 3.3224
-0.459 1.284 -0.654
22 22 22
ture.17,18 In this work, we have chosen the SPC/E model19 because it provides a good representation of the dielectric properties at ambient temperature as well as thermodynamic properties at a relatively low computational cost. The Na+ ion is modeled as a charged Lennard-Jones sphere20 using the parameters given in Table 1. The model and force field employed for the DS- anion are very similar to those employed in previous simulations of SDS monolayers,16,21 with some modifications taking into account the AMBER 94 force field.22 The hydrocarbon chain was modeled as a terminal united atom CH3 group connected to a hydrocarbon chain of 11 united atom CH2 sites. The five atoms in the hydrophilic headgroup SO4 were explicitly modeled. The bonds between the atoms were described with harmonic potentials U ) kb(r - req)2/2 with the parameters shown in Table 2. The bending and torsional forces in the DS- anion were described by the following potential energies:
Ubend )
kijk (θijk - θe)2 2
Utor ) Aijkl[1 + cos(mφijkl + δ)]
(1) (2)
where θijk is the angle determined by the atoms ijk, and φijkl is the dihedral angle between four consecutive atoms i, j, k, l. The constants m and δ were given by m ) 3 and δ ) 0.21 The values of the parameters in eqs 1 and 2 are also reported in Table 2. Each atom interacts with the atoms of other molecules via electrostatic and LennardJones interactions. The distribution of charges in the surfactants was chosen following previous work based on semiempirical quantum mechanics calculations.23 The group charges and the parameters for the Lennard-Jones interactions are shown in Table 1. These values are very similar to those employed in previous simulations.16,21 The only significant difference is the larger Lennard-Jones diameter of the S atom, chosen following ref 22. The simulations were carried out with the DLPOLY package24 running in the HPCx supercomputing center at Edinburgh, using 32-128 processors depending on the system size. All of the simulations were performed at ambient temperature (T ) 298 K) in the canonical (17) Guissani, Y.; Guillot, B. J. Chem. Phys. 1993, 98, 8221. (18) Bresme, F. J. Chem. Phys. 2001, 115, 7564. (19) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (20) Dang, L. X. J. Am. Chem. Soc. 1995, 117, 6954. (21) Schweighofer, K. J.; Essman, U.; Berkowitz, M. J. Phys. Chem. 1997, 101, 3793. Dominguez, H.; Berkowitz, M. L. J. Phys. Chem. 2000, 104, 5302. (22) Cornell, W. D.; et al. J. Am. Chem. Soc. 1995, 117, 7, 5179. (23) Shelley, J.; Watanabe, K.; Klein, M. L. Int. J. Quantum Chem: Quantum Biol. Symp. 1990, 17, 103. (24) Forester, T. R.; Smith, W. DLPOLY Package of Molecular Simulations v 2.13; CCLRC: Daresbury Lab, 2001.
Computer Simulation Studies of Newton Black Films
Langmuir, Vol. 20, No. 12, 2004 5129
Table 2. Bonding Interactions (Energies in kJ/mol and Distances in Å) interacting groups
interaction parameters
ref
bond CH3-CH2 bond CH2-CH2 bond CH2-O(ester) bond O(ester)-S bond S-O bending CH3-CH2-CH2 bending CH2-CH2-CH2 bending CH2-CH2-O(ester) bending CH2-O(ester)-S bending O(ester)-S-O bending O-S-O dihedral CH3-CH2-CH2-CH2 dihedral CH2-CH2-CH2-CH2 dihedral CH2-CH2-CH2-O(ester) dihedral CH2-CH2-O(ester)-S dihedral CH2-O(ester)-S-O
kb ) 2594.1, re ) 1.526 Å kb ) 2594.1, re ) 1.526 Å kb ) 2677.8, re ) 1.410 Å kb ) 2677.8, re ) 1.580 Å kb ) 7531.2, re ) 1.460 Å kijk ) 334.72, θe ) 109.5° kijk ) 334.72, θe ) 109.5° kijk ) 418.40, θe ) 109.5° kijk ) 418.40, θe ) 112.6° kijk ) 427.00, θe ) 102.6° kijk ) 427.00, θe ) 115.4° Aijkl ) 11.7152 Aijkl ) 11.7152 Aijkl ) 11.7152 Aijkl ) 9.6232 Aijkl ) 2.0920
22 22 22 22 21 22 22 22 22 22 22 21 21 21 21 21
ensemble using the Nose-Hoover thermostat,25 with a relaxation constant of 0.5 ps, and the equations of motion were solved with the Verlet leapfrog algorithm using a time step of 0.002 ps. The total simulation time was 1001 ps for each system, including 200 ps of equilibration. Our objective is to study the thin film properties as a function of film thickness. To obtain such films, we have performed simulations with different water contents, with the number of water molecules ranging from Nw ) 22NSDS (where NSDS is the number of SDS molecules) to Nw ) 0. In all of the simulations, we have considered 64 DS- anions per monolayer (NSDS ) 128). In our simulations, the films are surrounded by a vacuum in the direction normal to the bilayer plane, so that all of the simulations reported here correspond to zero normal pressure conditions (strictly speaking, the pressure in the exterior of the film is equal to the vapor pressure of water, which is negligible for all practical purposes). The dimensions of the simulation box in the bilayer plane were Lx ) Ly ) 46 Å, which corresponds to the experimentally determined surface area per surfactant4 33 Å2. Periodic boundary conditions were applied in all three directions. The cutoff for the van der Waals interactions was set to 14 Å, and the Coulombic interactions were computed using the Ewald summation method25 with a precision of 10-6-10-5. A version of the Ewald summation periodic only in the x and y directions was not employed because it considerably increases the CPU time and the simulation times become prohibitive (see, for instance, ref 26). At this point, it is important to note that in our simulations the periodic boundary conditions are also applied in the z direction, so the bilayer electrostatically interacts with its own images with a field proportional to the z component of the total dipole moment of the system (Pz).25,27 This interaction is present even if Lz is relatively large. Nonetheless, in our system, it is completely negligible because Pz ≈ 0 due to the symmetry of the system. To check the reliability of the selected Lz value, tests were performed with different configurations and different values of Lz to ensure that electrostatic interactions are not affected by our choice of Lz. Finally, we chose Lz between 100 and 120 Å depending on the simulation; in all cases, the selected value for Lz was much larger than the film size in the z axes. To minimize the equilibration time, the simulations were started from preequilibrated water slabs. First, we produced a film of equilibrated SPC/E water with ∼3000 water molecules at 298 K inside a box with Lx ) Ly ) 46 (25) Frenkel, D.; Smit, B. Understanding Molecular Simulations, 2nd ed.; Academic Press: London, U.K., 2002. (26) Spohr, E. J. Chem. Phys. 1997, 107, 6342. (27) Yeh, I.; Berkowitz, M. J. Chem. Phys. 1999, 111, 3155.
Figure 1. Evolution of the bilayer thickness, H(t), with simulation time with respect to the equilibrium thickness, Heq. Full line, Nw/NSDS ) 11.94, and dashed lines, Nw/NSDS ) 4.37.
Å and Lz ) 120 Å. Next, 64 DS- anions were added at each side of the water film, with the hydrophilic head in contact with the water interface, and 128 Na+ ions were randomly placed inside the water film. All of the water molecules overlapping with the newly inserted particles were removed. The resulting configuration was equilibrated until the film size and internal energy attained a stable value. The final configuration from this equilibration was employed as the starting configuration for the simulation involving Nw ) 22NSDS water molecules. Initial configurations for systems with smaller water content were obtained from the final configuration of the simulations of higher water content, by removing randomly selected water molecules. III. Results A. Film Thickness as a Function of Water Content. In this subsection, we discuss the film thickness, H, as a function of water content Nw/NSDS. During a typical simulation, the thickness decreases from the initial value (generated as described in the previous section) to a stable final value, H ) H(Nw/NSDS), assumed to be the equilibrium thickness of the film for the given number of water molecules. We have found that the relaxation of the bilayer from the initial configuration to the equilibrium thickness is rather slow. For equilibrium systems, we observe that the typical vertical displacement of a sulfate group after 1 ns is of the order of 〈∆z2〉1/2 ≈ 6.3 Å, which is rather small. Figure 1 illustrates the variation of the bilayer thickness with simulation time for two representative films. The thickness was defined as the average distance between the O(ester) atoms in surfactants located in
5130
Langmuir, Vol. 20, No. 12, 2004
Bresme and Faraudo
Figure 2. Dependence of film thickness, H, on water content. Nw/NSDS represents the number of water molecules per surfactant. The arrow indicates the location of the thikness, Hc; see text for details.
different monolayers. One can see that the relaxation to the equilibrium thickness involves times of the order of 200 ps. The slow motion toward the equilibrium thickness is a reflection of the slow atom dynamics in the NBF. We would like to emphasize that these simulations were started from configurations in which the water was already equilibrated. Earlier simulations of Newton black films were started from a lattice configuration, which is far from equilibrium, and the simulation time was restricted to 200 ps.16 According to the current investigation, these simulations were still too short to attain a proper equilibrium thickness. The equilibrium thickness was computed as an average over all of the configurations after the equilibration period. The dependence of the film thickness on water content is depicted in Figure 2. The solid line represents a simple theoretical prediction, assuming that the water inside the film has a constant density F ≈ 1 g/cm3:
H ) H0 +
mNSDS Nw FLxLyNa NSDS
(3)
where m/Na is the mass of one water molecule, F is the density of bulk water, 1 g/cm3, and Na is Avogadro’s number. We note that H0 is a free parameter independent of the water density. The slope of the theoretical line agrees with the simulations for film thicknesses H larger than Hc ≈ 7.5 Å (this corresponds to a water content of Nw/NSDS ≈ 2.25). Therefore, we expect that, for films with H > Hc, water inside the film behaves essentially as ordinary bulk water. It is interesting to note that, for thicknesses H < Hc, there is a change in the slope of the function H(Nw/ NSDS) (see Figure 2). At this point, the variation of the water thickness does not conform anymore to the idea of the slab having the water bulk density. We believe that this change in slope signals the region where the film completely loses the bulk water slab. The remaining water molecules then would be taking part in the solvation of the surfactant headgroups. Indeed, an inspection of representative snapshots for bilayers with different water content supports this view (cf., Figure 3a-c). One can barely see any free water in the middle of the film for very thin films (cf., Figure 3c). Even for thicker films (see Figure 3b), the distinction between the two monolayers becomes difficult in some regions. It is very interesting to note that the critical thickness, Hc, is close to the thickness of the NBF spontaneously formed in the experiments reported in ref 4.
Previous computational studies16 have reported a linear behavior in the dependence of the film thickness on water content down to zero water molecules per surfactant; that is, no change in the slope of the thickness was observed. This is in contrast with the results reported here. Analysis of Figure 9 in ref 16 shows that the simulation results are not compatible with the expected asymptotic result for the thickness of the film, as obtained by considering the bulk water model employed in their simulations. We believe this might indicate lack of equilibration in the simulations. As we have noted before, we find that very long simulations are needed to reach equilibration. One of the main strengths of computer simulations is the possibility of investigating in detail the structure of the black films. As mentioned in the Introduction, a detailed knowledge of the structure of the film is of major importance to understand the forces acting in the system. The snapshots presented in Figure 3 provide a qualitative idea of the structure of the system, but in the following subsections we report a more quantitative description. B. Density Profiles. The changes in the water structure with the film thickness can be investigated in more detail by computing the density profiles. Figure 4 reports the water density profile for different bilayers. The profiles corresponding to the left and right monolayer have been averaged. We would like to point out that in our simulations the density profiles are not in general exactly symmetric with respect to the middle plane (z ) 0). For a very long simulation, we would expect a symmetrical density profile, but our simulations illustrate the notion that some asymmetries can survive during nonnegligible time intervals, at least 1 ns in our case. The water density profiles (see Figure 4) clearly show a film with a well-defined bulk region for the larger systems with a broad interfacial region. For example, for Nw/NSDS ) 11.94, we estimate that the thickness of the watersurfactant monolayer interface is ∼10 Å. As a comparison, the width of the water-air interface at ambient temperature, as obtained from experiments, is 4 Å.28 An increase in the water-surfactant interfacial thickness is expected because the surface tension of this interface is smaller than that in the water-air case. As an example of a system with an intermediate water content, let us consider the case Nw/NSDS ) 4.37 (see Figure 3b). In this bilayer, the lack of a clear water slab is noticeable, even though the profile suggests there is still some water in the middle of the NBF with a density which is very close to that of bulk water. For even thinner bilayers, Nw/NSDS ) 2.14, almost all of the water in the film is part of the bilayer structure, and clearly there is no bulk water in the middle of the film (cf., Figure 3c). These results indicate that the film loses all of the bulk water at a specific thickness, and the remaining water should be part of the bilayer structure. These observations are compatible with the results presented in Figure 2 and help one to understand the distinctive change in slope observed in the dependence of the film thickness on water content. Figure 5 represents the density profiles for the CH3, the sulfur atom in the surfactant headgroup, as well as the Na+ ions and the water oxygen. The profiles indicate that the reduction in the bilayer thickness does not result in qualitative changes in the density profiles, except those due to the overlapping of density distributions corresponding to the two monolayers. The density profiles also show that water penetrates the hydrocarbon side, solvating the surfactant headgroups. Indeed, it is very clear that there is a change in the shape of the water-oxygen (28) Matsumoto, M.; Kataoka, Y. J. Chem. Phys. 1988, 88, 3233.
Computer Simulation Studies of Newton Black Films
Langmuir, Vol. 20, No. 12, 2004 5131
Figure 3. Representative snapshots of thin films with different water contents: (a) Nw/NSDS ) 11.94, (b) Nw/NSDS ) 4.37, and (c) Nw/NSDS ) 2.14.
Figure 4. Oxygen(water) density profiles across the bilayer film as a function of film water content.
density profile appearing around the maximum of the sulfur profile as a distinctive shoulder. The number of water molecules nwh “inside” the monolayer structure in the “hydrocarbon” side of the monolayer can be estimated from:
nwh ) LxLy
∫zz
max
min
Fwater(z) dz
(4)
where zmin represents an arbitrary position in the hydro-
carbon region where water molecules do not penetrate, and zmax represents the position at which the sulfur density profile is maximum. The results computed from eq 4 for the systems with Nw/NSDS ) 11.4, 4.37, and 2.14 are nwh ) 72, 68, and 80 molecules, respectively. We expect that the surfactants are similarly solvated on the opposite side (the “water side”) of the monolayer, but this is not easy to evaluate from the density profiles because there are contributions coming from the water bulk molecules. For large monolayer-monolayer separations, H > Hc, we assume that the number of water molecules coating the “water side” of the monolayer is the same as that in the hydrocarbon side. In this case, there would be around 2 × 70 water molecules per monolayer or, in other words, 2.2 water molecules per surfactant. This shows that there is a considerable number of water molecules taking part in the monolayer structure. Also, it is very interesting to realize that the change in slope in the function H ) H(Nw/ NSDS) appears at Nw/NSDS ≈ 2.25, which is consistent with this estimate. Hence, one can conclude that the critical value Hc ≈ 7.5 Å corresponds to a system with a water content consisting only of solvation water. Our interpretation is that the water molecules in excess from the quantity Nw/NSDS ≈ 2.25 behave as bulk water and would not take part in the monolayer structure. Finally, let us discuss the distribution of the Na+ ions, also reported in Figure 5a-c. The vast majority of Na+
5132
Langmuir, Vol. 20, No. 12, 2004
Bresme and Faraudo
Figure 5. Density profiles across the bilayer film for CH3 and sulfur groups, oxygen(water), and sodium. (a) Nw/NSDS ) 11.94, (b) Nw/NSDS ) 4.37, and (c) Nw/NSDS ) 2.14.
ions are close to the DS- monolayers. The profiles show there is an important overlap between the distributions of Na+ ions and S atoms. This structure departs from that considered in different works,8,29 where the counterions and the surfactant heads are in different planes. These density plots show that the distribution of the S atoms is quite broad as compared with the size of the film, so the adsorbed layers cannot be modeled as a simple plane. Also, our results indicate that most of the ions are taking part in the structure of the layer. This implies that an important fraction of the Na+ ions are very likely to be adsorbed in the monolayer, or, in other words, the surfactants are not fully ionized, on average. A rough estimate shows that nearly 90% of the Na+ ions are adsorbed in the case of films with H > Hc. This result should have important implications on the stability of the film, because the distribution of ions determines both the electrostatic interactions inside the film and the osmotic pressure. A detailed analysis of the electrostatics for the present model is in progress, and it will be reported elsewhere.30 C. Roughness. The roughness of a fluid-fluid interface is an important structural quantity, reflecting a compromise between the elastic properties of the interface and thermal fluctuations in the form of capillary waves. It can be computed from the height-height correlation function, or roughness function,
ζ(|R B |) ) x〈(z(r b) - z(r b+R B ))2 〉
(5)
where z(r b) and z(r b+ R B ) are the z coordinates of two S atoms in the same monolayer. For systems above the socalled roughening temperature, the quantity ζ diverges
Figure 6. Monolayer roughness for Nw/NSDS ) 11.94 bilayers with 128 and 512 surfactant molecules (full lines). The dashed line represents the system Nw/NSDS ) 4.31, and the dotteddashed line represents the roughness of the Nw/NSDS ) 11.94 system using different diameters for sulfur and sodium. See text for details.
for R f ∞ as ζ(R) ∝ ln(R), so that the film is destabilized by the capillary waves induced by thermal fluctuations. For stable films, the roughness ζ(R) approaches a constant value as R f ∞. The roughness is reported in Figure 6 for two films with water content Nw/NSDS ) 11.94 and 4.37. It is interesting to note that the roughness for these two systems is essentially the same. To compute the asymptotic value of the roughness, we have performed additional simulations using a larger simulation box Lx ) Ly ) 96 Å containing 256 surfactants per monolayer and Nw/NSDS
Computer Simulation Studies of Newton Black Films
Langmuir, Vol. 20, No. 12, 2004 5133
Table 3. Average Hydrocarbon Chain Length and Tilt Angle for Different Thin Bilayer Films Nw/NSDS
lcc (Å)
cos(θ)
11.94 4.37 2.14
10.5 ( 0.1 10.5 ( 0.1 10.6 ( 0.1
0.79 ( 0.01 0.80 ( 0.01 0.80 ( 0.01
) 11.94 water molecules per surfactant. These large-scale simulations show that the roughness function reaches an asymptotic value ζ(R f ∞) ≈ 2.5 Å, which is in very good agreement with the experimental estimate, 2.7 ( 0.1 Å,4 and it gives further support to the model used in this work. Theoretical calculations, however, predict a roughness of 2.4 Å,31 which is slightly smaller. We note that these calculations assumed a bending rigidity of kBT and a monolayer surface tension of ∼76 mN/m, that is, that of pure water, which is twice as large as the surface tension of the corresponding monolayer. At this point, it is interesting to check the sensitivity of our results to the parameters employed in the model, especially to the diameters of the sulfur and sodium atoms. We have performed an additional simulation using different values for the Lennard-Jones diameters of these atoms, employing the values used in ref 21: σS ) 3.55 Å and σNa ) 2.275 Å. All of the other parameters are the same as those considered before. The new results (cf., Figure 6) show that the roughness is very sensitive to small changes in the interatomic potentials. Our results indicate that a reduction of the cation and sulfur diameters results in an increased roughness, up to 20%, overestimating the experimental result. Therefore, it seems that the roughness is a rather sensitive property that could be useful in the assessment of the accuracy of interatomic potentials. D. Tilt Angle. Table 3 contains results for the average hydrocarbon chain length and tilt angle for three different film thicknesses. The average chain length was computed from the average distance between the terminal CH3 and the CH2 group bonded to the O(ester) atom, and the tilt angle from the average angle between the vector connecting these two groups and the axes normal to the bilayer plane. The results are reported in Table 3. The chain length and tilt angle are rather insensitive to the film thickness. This is compatible with the experimental observation that the aliphatic structure in the NBF is the same for different film thicknesses.4 This result justifies the approximation taken in several theoretical treatments,8,9 where the hydrophobic region of the film is decoupled from the aqueous core. The hydrocarbon chain lengths obtained from our simulations are in agreement with the estimates obtained for SDS monolayers at the water interface.21 The tilt angle in the present model is 37°, which is slightly smaller than that reported in other simulations, approximately 40° in ref 32 and 45° in ref 16. Hence, in our case, the hydrocarbon chains are in a more straight conformation. Other simulation results for SDS monolayers at the water/vapor interface indicate that the distribution of tilt angles is very broad, and no clear tendency is observed for the tilt angle of the chains.21 E. Electron Density Profile. To compare our results with X-ray experimental data, we show in Figure 7 the electron density profile for the systems with Nw/NSDS ) 11.94, 4.37, and 2.14. The basic difference between these (29) Gulbrand, L.; Jo¨nsson, B.; Wennersto¨m, H.; Linse, P. J. Chem. Phys. 1984, 80, 2221. (30) Faraudo, J.; Bresme, F. Phys. Rev. Lett. 2004, in press. (31) Daillant, J.; Bosio, L.; Benattar, J. J.; Meunier, J. Europhys. Lett. 1989, 8, 453. (32) Bareman, J. P.; Klein, M. L. J. Phys. Chem. 1990, 94, 5202.
Figure 7. Electron density profile for the SDS bilayer for systems with different water contents. The experimental electron density profile reported in ref 5 has also been depicted for the sake of comparison.
results and those obtained in the simulations of ref 16 is that we do not obtain the unexpected peak in the middle of the profile, which was probably a signal that the aqueous core of the simulated film was not in equilibrium. The central region of the film in the case of the Nw/NSDS ) 4.37 system agrees fairly well with the experimental result of ref 4. However, our results predict that the hydrocarbon chains do not show a well-defined bulk region. This discrepancy with the experiments has also been observed in previous computer simulations.16 We note that our results are consistent with the fact that the hydrocarbon chains are rather short and they undergo thermal fluctuations along the z axis. The origin of the differences observed between the experimental and simulation results remains to be clarified. F. Pair Correlation Functions. The pair correlation functions provide a more detailed image of the structure of the atoms in the film. Figure 8a,b shows the functions FRβ(r) for S-S and S-Na and Na-Na correlations in films with water content Nw/NSDS ) 11.94, 4.37, and 2.14. FRβ(r) is the average density of atoms of species, β, lying a distance r from an atom of species R. We represent this function instead of the normalized radial distribution function, g(r), for the sake of clarity. The correct normalization of g(r) for these systems (see, for instance, the definition of g(r) in ref 33), in terms of the density of the corresponding ideal gas, depends on the box volume chosen in the simulation. To avoid this ambiguity in the representation of g(r), we have considered the density FRβ which is independent of the volume of the simulation box. In Figure 8, we represent FRβ for S atoms. It shows that there is a considerable ordering of the headgroups in the film plane. This is particularly clear in the thicker film (Nw/NSDS ) 11.94), where there is no possibility of monolayer overlap and the structure is due only to headgroups in the same monolayer. We find that the structure is enhanced as the water content decreases, particularly the second peak. The average thickness of the film for the case Nw/NSDS ) 4.37 is H ≈ 11 Å, so this increase is likely due to correlations between atoms in different monolayers. Figure 9a gives an alternative view of the monolayer structure for the system Nw/NSDS ) 11.94. Interestingly, the headgroups are not distributed uniformly over the surface. Rather, we find that the monolayer consists of large headgroup clusters and large cavities filled with water. This image clearly departs from a model (33) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1993.
5134
Langmuir, Vol. 20, No. 12, 2004
Bresme and Faraudo
Figure 8. Pair correlation functions, FRβ (see text for definition), for bilayers with different water contents. (a) Sulfur-sulfur, (b) sodium-sodium, (c) sulfur-sodium.
Figure 9. Snapshots of representative monolayers for films with (a) Nw/NSDS ) 11.94 and (b) Nw/NSDS ) 4.37.
in which the surfactants are uniformly distributed in the monolayer and suggests that the systems may undergo relatively large fluctuation in the surfactant density. We would like to point out that the existence of cavities in ionic systems has been observed and characterized very recently by one of us in the case of ionic liquids close to the coexistence region.34 We observe the formation of
cavities also in films with different water content, see, for instance, the system with Nw/NSDS ) 4.37 (cf., Figure 9b). As far as we know, this observation has not been reported before. The consequences that this hole nucleation may have on the electric field patterns generated at short (34) Bresme, F.; Alejandre, J. J. Chem. Phys. 2003, 118, 4134.
Computer Simulation Studies of Newton Black Films
Langmuir, Vol. 20, No. 12, 2004 5135
Figure 10. (a) Lateral diffusion coefficient as a function of water content. The dots correspond to the simulation results, and the line corresponds to eq 7. (b) Mean squared displacement in the direction normal to the bilayer plane as a function of time for different systems. The solid, dashed, and dotted lines correspond to Nw/NSDS ) 11.94, 4.37, and 2.14, respectively.
interbilayer distances and the film stability require further investigation. The function FRβ for Na atoms (cf., Figure 8b) also shows distinctive features and illustrates much more clearly the structural changes induced by a decrease in the bilayer thickness. The first peak indicates that the sodium ions undergo a reorganization as the water content decreases. An interpretation of the structure attained by the film for the systems with smaller water content is not trivial because there is a large overlap between both layers, but, on the whole, looking at the S-S and Na-Na correlation functions, the Na-Na one shows the largest changes, indicating that the counterions undergo a larger reorganization than do the headgroups. Our results for the S-S correlation functions are rather different from the pair correlation functions reported by Gamba and co-workers.16 Some differences between both sets of results are expected because our model is slightly different. Nonetheless, for a monolayer in which the headgroups show a large degree of in-plane ordering, we think that the S-S correlation function should exhibit a marked structure, like the one we obtain here. We note that our results for S-Na (cf., Figure 8c) and Na-Na (cf., Figure 8b) show similarities to the results reported by previous authors.16 For instance, S-Na exhibits the characteristic splitting of the first peak at 3.8 Å due to counterions located behind the headgroup and on the terminal O(ester) atoms. The S-Na coordination number for NBF with water content Nw/NSDS ) 11.94, 4.37, and 2.14 is 2.6, 2.7, and 3.4, respectively. The significant increase in the coordination number is due to the overlap of the two monolayers when the NBF water content is very small. As a test to address the sensitivity of the correlation functions to the interatomic potentials, we have performed
a simulation in which we have changed the diameter of the sulfur and sodium atoms. We have chosen as new Lennard-Jones diameters the values employed in ref 21 (σS ) 3.55 Å and σNa ) 2.275 Å), maintaining all of the other details of the simulation. The results from this modified model are presented in Figure 8a-c. We essentially observe the same structure for the S-S correlation, which appears to be much less sensitive to the changes than are the Na-Na and S-Na pair correlation functions. The new bilayer presents qualitatively the same structure, although the fact that now the sodium counterion is much smaller results in a decrease in the S-Na and Na-Na distances. G. Water Mobility. In the following, we analyze the water diffusion inside the thin film. This dynamic information helps one to understand the role played by water in the film and complements the more static picture given by the analysis of density profiles and other structural properties. Due to the anisotropy of the system, one has to distinguish between the lateral motion of the molecules in the xy plane and the vertical motion in the z direction. In the xy plane, ordinary diffusion is observed with a substantial reduction in the diffusion coefficient as compared with that of bulk water. The mean squared displacement in the XY increases linearly with time, as expected in ordinary Fickian diffusion. The lateral diffusion coefficient D || is computed from our simulations by using the relation:
〈(x(t) - x0)2 + (y(t) - y0)2〉 ) 4D|t
(6)
The diffusive behavior described by this relation is reached in our simulation for times larger than 1 ps. The evolution of D| with the water content of the film is shown in Figure
5136
Langmuir, Vol. 20, No. 12, 2004
Bresme and Faraudo
10a. It is clear that the lateral mobility of water molecules is highly reduced by the adsorbed surfactant layers even in the case of relatively high water contents (corresponding to larger separations between both monolayers). This reduction has different sources. The first one is the fact that some water molecules are immobilized because they are solvating the adsorbed layers (i.e., they are inmobilized by the bilayer structure). This effect can be estimated easily by assuming that the reduction in the diffusion coefficient is proportional to the fraction of these immobilized water molecules:
[
Dth | ) D0 1 -
]
4nwh Nw
(7)
In eq 7, the number of inmobilized water molecules is assumed to be approximately equal to 4 times the number of water molecules nwh which are solvating the surfactant layer in the hydrocarbon side (see eq 4), and D0 = 0.25 A2/ps is the diffusion coefficient of bulk SPC/E water. The estimate from eq 7 is also shown in Figure 10a (solid line). It is clear that water inmobilization cannot explain the observed reduction in water mobility. Another factor that could reduce the lateral diffusion of water molecules is the headgroup protrusion inside the film. This could make difficult the mobility of the water molecules; however, it is difficult to make an estimation of this effect in our system. Another factor that could influence the mobility is the water structuring within the bilayer core. An analysis of this structure is currently under investigation. The mobility of water in the z direction presents a complex behavior. Figure 10b shows the displacement 〈(z(t) - z0)2〉 for films with different water content. The mean squared displacement increases with time up to a maximum value, which depends on the bilayer thickness. This behavior in the mean squared displacement can be rationalized considering a particle diffusing within a confined region of constant thickness L, that is, within a pore. For sufficiently large times, the particle explores the whole system and the mean squared displacement reaches a limiting value 〈(z(t f ∞) - z0)2〉 f L2/12. In the case of our bilayer, the situation is more complex, because there are water molecules which are part of the monolayer structure and therefore exhibit a different dynamics than the water molecules in the water slab. Also, the shape of the bilayer fluctuates, so the concept of a well-defined volume where the molecules can diffuse does not strictly apply to our system. In any case, the results suggest that the water molecules diffuse in the vertical direction as if they were confined to a region smaller than that predicted by L2/12, where L can be estimated from the water density profile as the distance between the two points where the profile becomes nearly zero. For instance, for Nw/NSDS ) 4.37, we obtain 〈(z(t) - z0)2〉 ) 24 Å2, as compared with L2/12 ) 292/12 ) 70 Å2 (see Figure 5b). It is remarkable that the maximum vertical displacement of water molecules is very small, only a few water molecule diameters. IV. Conclusions and Final Remarks In this work, we have studied via molecular dynamics simulations the structure and dynamics of a sodium dodecyl sulfate Newton black film, with a surface area of 33 Å2 per surfactant, corresponding to the surface area found in X-ray experimental studies of SDS black films.4 One of the conclusions from our work is that the simulations of thin bilayer films are computationally demanding, and simulation times of the order of the nanosecond are needed to attain proper mechanical equilibrium states.
Computer simulations of SDS black films with different water contents at constant surface area reveal that the film thickness varies linearly with water content for H > Hc ) 7.5 Å and indicate that the density of water inside the film essentially is that of bulk water. We find that the thickness exhibits a distinctive change in slope for a bilayer with water content of 2.25 water molecules per surfactant. We conclude from the analysis of the density profiles that this discontinuity is related to the loss of bulk water in the film. Indeed, we estimate from our simulations that the amount of water taking part in the bilayer structure as solvation water is ∼2.2 water molecules per surfactant. Interestingly, this value of 2.2 water molecules per surfactant corresponds to an average O(ester)-O(ester) distance of ∼7.6 Å. The thickness obtained for a freestanding NBF at atmospheric conditions in the experiments4 is ∼7.5 Å. Comparison of our X-ray electron density profiles with the experimental one shows that the bilayer that better fits the experimental profile in the central region should correspond to a composition between 2.14 and 4.37 water molecules per surfactant. These values are slightly lower than those predicted in previous simulations,16 4-6 water molecules per surfactant. A proper prediction of the equilibrium thickness observed in these experiments requires an accurate calculation of the free energy of the bilayer as a function of the thickness. Because the aliphatic region is rather insensitive to the bilayer thickness, we expect that the free energy associated with the hydrocarbon chains should play a minor role in determining the equilibrium thickness of the NBF. We believe that the main contributions should come from the ionic interactions and the reorganization of the water structure. We have also computed the roughness of the monolayers for bilayers with relatively large and medium water contents, 11.4 and 4.37 water molecules per surfactant. Our model predicts a roughness of 2.5 Å which compares favorably with the experimental and theoretical estimates, 2.7 and 2.4 Å, respectively. We have found that the roughness is very sensitive to the interatomic potentials. We have shown that a small change in the sulfur and sodium Lennard-Jones diameters results in large changes in the roughness, in our case up to a 20% increase. Therefore, the roughness reveals itself as a property that could be used to assess the accuracy of semiempirical potentials. Interestingly, we have found that the monolayers in the NBF show a high degree of in-plane ordering; this is revealed by the sulfur-sulfur pair correlation functions and representative snapshots from the simulations. We find that the monolayer consists of clusters of headgroups and cavities in which only water is present. This image departs from that of a homogeneous distribution of surfactants in the monolayer plane. These structures could play an important role in, for instance, determining the interactions between monolayers at short distances, where electrostatic forces are expected to become increasingly important and the discreteness of the monolayer structure may be relevant in determining the final bilayer electric field. Also, the nucleation of holes in the NBF should be relevant to the film stability. Clearly, these observations and their implications in the physics of bilayer films require further investigation. Finally, we have analyzed the water mobility in the film as a function of water content and orientation, that is, parallel and perpendicular diffusion. In general, we find that the dependence of the diffusion coefficient on water content does not have a trivial explanation. Our results indicate that the slowing down of the water mobility
Computer Simulation Studies of Newton Black Films
at intermediate/small thickness cannot be explained only in terms of the expected amount of free water at a specific thickness. The mean squared displacement of water in the direction normal to the bilayer plane is more complex than that in the lateral direction. Qualitatively, this behavior can be explained considering the movement of water molecules within a pore consisting of two parallel walls. Our results indicate that the movement in the vertical direction is confined to a very small region, much smaller that the actual bilayer thickness measured as the O(ester)O(ester) distance. According to our data, the expected
Langmuir, Vol. 20, No. 12, 2004 5137
lateral diffusion of water at the equilibrium film thickness (between 2.14 and 4.37 water molecules/surfactant) should be very small, between 0.01 and 0.04 Å2/ps; that is, the water molecules should be essentially frozen in the equilibrium state. Acknowledgment. This work has been supported by EPSRC Research Grant No. GR/R39726/01. Computer time in the HPCx through the Materials Chemistry Consortium (U.K.) is gratefully acknowledged. LA036026W
