Liquid and

5302

J. Phys. Chem. B 2000, 104, 5302-5308

Computer Simulations of Sodium Dodecyl Sulfate at Liquid/Liquid and Liquid/Vapor Interfaces H. Dominguez and M. L. Berkowitz* Department of Chemistry, UniVersity of North Carolina at Chapel Hill, Chapel Hill, North Carolina 27599 ReceiVed: December 23, 1999; In Final Form: March 15, 2000

We performed molecular dynamics simulations on sodium dodecyl sulfate (SDS) monolayer at the water/ carbon tetrachloride and water/vapor interfaces. We observed that the tails are slightly less ordered when SDS is at the water/vapor interface. Also, at the water/carbon tetrachloride interface the amphiphilic molecule is less inclined with respect to the surface normal compared to the one at the water/vapor interface. We also carried out investigations of the electrostatic properties and surface tension of the water/carbon tetrachloride interface in the presence of the SDS monolayer with and without the dipole correction term in the Ewald sum. We observed that an inclusion of the correction term makes a difference in quantitative results.

1. Introduction The behavior of amphiphilic molecules at liquid/vapor and liquid/liquid interfaces has been a matter of scientific interest for a long time. A large variety of modern experimental techniques including fluorescence, resonance Raman scattering, neutron reflection, second harmonic generation, vibrational sumfrequency spectroscopy, time-resolved quasi-elastic laser scattering are used to study structural and dynamical properties of these systems.1-10 At the same time, due to the substantial increase in the computational power, computer simulations became an important tool for the study of complex interfacial systems on a detailed molecular level (e.g., refs 11-20). However, the results of the simulations may be very sensitive to the details of the simulations, i.e., the potential model (force field) used, the choice of an appropriate boundary conditions and geometry, accurate treatment of long-range forces, etc. Recently we performed a molecular dynamics study of one of the most common surfactantssthe anionic detergent sodium dodecyl sulfate (SDS)sat the water/vapor and at the water/ carbon tetrachloride interface.13,14 We found that at low surface coverage the SDS molecular configuration at the water/vapor interface is very different from the one at water/carbon tetrachloride interface.13 At the water/vapor interface the surfactant molecule is bent on average, giving rise to two domains within the molecule. The first domain containing the head group and several methylene groups was solvated in water, while the second domain containing the rest of the molecule lay down on the water surface. In contrast to this result the SDS molecule at the water/carbon tetrachloride interface was straight on average, with an inclination of approximately 40° from the interface normal. Our next study of SDS monolayer at the water/ carbon tetrachloride interface (at a coverage of 45 Å2/molecule) was focused on the structural and dynamical properties of water and on the electrical properties such as dipole potential of the interface.14 In this work we continue our study of SDS monolayer at a surface coverage of 45 Å2/molecule and compare the average structure of a surfactant molecule at water/vapor and water/carbon tetrachloride interfaces. (The value of 45 Å2/ molecule was chosen because it is equal to the area at the critical * Corresponding author. E-mail: [email protected].

micelle concentration for the SDS molecule at the water/vapor interface, as found in the neutron reflection experiment.2 At the CCl4/water interface a saturation area of 59 Å2/molecule was reported in the literature5). We also revisit the issue of the value of the dipole potential in view of the recent work from our laboratory on the treatment of long-ranged electrostatic forces for systems in slab geometry.21 Finally, we calculate the change in the value of surface tension when SDS is added to the interface. 2. Computational Method and Model For the present study we performed six different simulations on systems that contained different interfaces. Two of them were performed on an SDS monolayer at the CCl4/water interface. We used the SPC water model22 in the simulations. Most of the parameters for the carbon tetrachloride potential were the same as the ones used in our previous simulations, except that now we treated the molecules as rigid. Bond lengths were constrained using the SHAKE algorithm with a tolerance of 10-4. The amphiphilic molecule was considered as a hydrocarbon chain of 12 united carbon atoms attached to a head group SO4. The head group atoms were explicitly modeled. Except for the dihedral potential in the head group, the parameters for the SDS were the same as used in our previous work.13 We changed the functional form of the head group dihedral potential from the Ryckaert-Bellemans to the one used in AMBER with parameters from AMBER,23 since we think the latter better describes the torsional motion of the head groups. The parameters for the head group torsional potential used in this work are given in Table 1. All simulations were carried out in the NVT ensemble with a time step of 0.002 ps using the DL-POLY package.24 The temperature was controlled using the HooverNose thermostat with relaxation time of 0.2 ps.25 All simulations were performed at T ) 300 K. For the long-range electrostatic potential we used the particle mesh Ewald method with a precision of 10-4. The van der Waals interactions were cut off at 10 Å. The procedure for the preparation of the initial configuration was performed similarly to the one described in ref 13. The monolayer containing 36 SDS molecules in all-trans configuration was placed into a rectangular box with x and y dimensions

10.1021/jp994479x CCC: $19.00 © 2000 American Chemical Society Published on Web 05/13/2000

SDS at Liquid/Liquid and Liquid/Vapor Interfaces

J. Phys. Chem. B, Vol. 104, No. 22, 2000 5303

TABLE 1: Parameters for the SDS Head Group Dihedral Potential group

torsion (kcal mol-1)

γ(rad)

n

C-C-C-O C-C-O-S C-O-S-O

1.000 0.725 0.250

0 0 0

3 3 3

of 40.249 Å, corresponding to an area of 45 Å2 per head group. The z dimension of the box was set to 150 Å. The box is longer in the z direction to accommodate two liquid slabs and to prevent the formation of a second water/CCl4 interface. Instead, two vapor/liquid interfaces at the opposite ends of the box (on the z < 0 side of the vapor/water interface and on the z > 0 side of the liquid CCl4/vapor interface) were present in the box. The first simulation of the water/SDS/CCl4 system was performed without the correction for the total dipole moment of the system (this correction term which should be included in simulations with the slab geometry is described in refs 21 and 26; it is also briefly described below); the second simulation included the dipole correction term. Initially we performed a short MD simulation on SDS molecules at T ) 300 K with the head groups of the molecules pinned; then we increased the temperature to T ) 400 K and ran the system for 20 ps in order to randomize the tails. Subsequently we decreased the temperature by performing several short runs, each one for 20 ps, until we reached T ) 300 K. After this equilibration we added 1185 water molecules surrounding head groups and a layer of 415 CCl4 molecules was placed in the region of the tails. Water was placed at z < 0 and CCl4 at z > 0. The system was then equilibrated for 100 ps. Finally, 36 sodium anions were randomly inserted in the interfacial region. This final configuration was then equilibrated for 550 ps, and the data for analysis were collected from a 1 ns run. Configurational energy was monitored as a function of time to make sure that the system reached equilibrium. We also monitored the distribution of molecules in the system: particularly we examined the number of CCl4 molecules in the region of SDS tails and observed that this number did not change considerably during the simulations. A third simulation was performed on an SDS monolayer at the water/vapor interface in order to compare the results with the ones from the simulation on the monolayer at the CCl4/ water interface. This system was constructed from the SDS monolayer at the CCl4/water interface by removing CCl4 molecules. Again the system was equilibrated for 500 ps and then we ran the system for another 1 ns for data analysis. We also performed a simulation on the vapor/water/CCl4/ vapor system in order to see how the width of the water/CCl4 interface changes when the surfactant is added to the system. Finally simulations of vapor/water/vapor and vapor/carbon tetrachloride/vapor systems were also performed in order to get values for the surface tension at water/vapor and CCl4/vapor interfaces. All these simulations were performed for a duration of 1 ns. 3. Results 3.1. SDS Configuration, Head Group, and Hydrocarbon Chain Locations. In our previously published work on SDS monolayer at the water/vapor and water/carbon tetrachloride interfaces, we concentrated on the study of water structure. In this work we want to focus on the study of SDS and also address the issues that were not discussed in our previous work. In Figure 1 we show the z-dependent density profiles for the liquids and head groups and hydrocarbon tails of solute molecules. The

Figure 1. Density profile for the of SDS monolayer at the vapor/ water interface (top) and at the CCl4/water (bottom). Water is depicted by the solid line, CCl4 by the dotted line, SDS head groups by the dashed line, and SDS tails by the dash-dotted line.

head group density profiles include the SO4- group and the Na+ counterion. Density profiles shown in Figure 1 are very similar to those obtained in our previous work.14 To get more information on the chain configurations of the SDS molecules, we measure their total chain lengths, which are defined as the distance, d, from the first carbon to the last carbon of the chain, resulting in a value of 11.6 ( 0.1 Å. A rough estimate of the chain tilt can be obtained from the projection of this distance along the normal to the surface, dz (cos θ ) dz/d). For the monolayer at water/CCl4 interface we find that the tilt angle is 39° ( 3° which is similar to the one calculated in a previous simulation of the same system at a much lower surface coverage.13 For the case of the monolayer at the water/vapor interface the values for the length of the chain and the tilt angle are 11.4 ( 0.1 Å and 47° ( 2° respectively. These results suggest that the tails at the CCl4/water interface are in a somewhat more straight configuration than at the water/vapor interface, as was observed in our previous simulations.13 More detailed information on the positions of the carbon atoms of the chains can be obtained from probability distributions of each carbon. These results are shown in Figure 2. The leftmost distribution is for the oxygen attached to the first CH2 group in the tail, and the rightmost distribution is for the last CH3 group. For the monolayers at CCl4/water interface, the average position of each site seems to change smoothly without

5304 J. Phys. Chem. B, Vol. 104, No. 22, 2000

Dominguez and Berkowitz

Figure 3. Cosine of the angle between the C1-Cn (n ) 2, 3, ..., 12) vector and the vector normal to the interface. Circles are for the SDS at the CCl4/water interface, and squares are for SDS at the water/vapor interface.

Figure 2. z-Dependent probability distribution for each carbon in the chain and for the oxygen attached to the tail. The top picture is for the SDS at the water/vapor interface, and the middle picture is for the SDS at the CCl4/water interface. The leftmost distribution refers to the oxygen attached to the first carbon (C1) in the tail, and the rightmost distribution refers to the last carbon in the tail (C12). The bottom picture shows the average z-position of each carbon distribution from the top and the middle panels. Circles are for the SDS at the CCl4/water interface, and squares are for the SDS at the water/vapor interface.

any sharp variation. This result suggests that the carbons are regularly spaced along the surface normal. For the monolayer at the water/vapor interface we also observe smooth changes, but the spacing between the positions of the peaks in the distributions is smaller compared to the ones at the CCl4/water interface. This again indicates that the surfactants have a larger inclination at the water/vapor interface. Finally, on Figure 3 we display the average values of the angle between the C1-Cn vector (n ) 2, 3, ..., 12) and the normal to the interface. From this figure we notice that, for the monolayer at CCl4/water interface, values for this angle for carbons down the chain (starting at the fifth carbon) reach a plateau, suggesting that those carbons have the same inclination. This inclination (≈39°) is in good agreement with our previous observation. No plateau is observed for the monolayer at the water/vapor interface, which indicates that chains are more disordered in this case. The angle distribution shown in Figure 3 is consistent with the previously

Figure 4. SCD order parameter as a function of the carbon position. Circles are for the SDS at the CCl4/water interface, and squares are for SDS at the water/vapor interface.

stated conclusion that the SDS chains are more tilted at the water/vapor interface. 3.2. Order Parameters and Number of Gauche Defects. The ordering of the tails in phospholipid membranes is usually characterized by the so-called deuterium order parameter, SCD, which shows the average inclination of the C-D bond with respect to the bilayer normal. To measure these, the hydrogens of the chains are selectively replaced by deuteriums and the NMR technique is used. In computer simulations that use united CHn atoms the SCD order parameter is calculated using the following formula:27

SCD(2/3)Sxx + (1/3)Syy

(1)

Sij ) (1/2)〈3 cos θi cos θj - δij〉

(2)

where

i, j ) x, y, z, and θi is the angle between the ith molecular axis and the normal to the interface (see the details in ref 27). In Figure 4 we show the results obtained for the SDS monolayers at both interfaces. For the SDS at the CCl4/water interface we observe a plateau region extending through carbons 3-6. The SCD parameter is decreasing starting from carbon 7,

SDS at Liquid/Liquid and Liquid/Vapor Interfaces

J. Phys. Chem. B, Vol. 104, No. 22, 2000 5305

Figure 5. Probability of gauche defects as a function of carbon position. Circles are for the SDS at the CCl4/water interface, and squares are for SDS at the water/vapor interface.

indicating that last carbon atoms of the chain are distributed in a more isotropical way. For the water/vapor interface no plateau is observed and the order is always lower when compared to the the CCl4/water interface. It is also possible to characterize the distribution of chain conformations by the probability of gauche defects (Figure 5). From this figure we observe that the S-O-C1-C2 dihedral is always almost trans, while the following dihedral angle has only 60% trans character. Note that this value is nearly the same for simulations at the water/vapor and CCl4/water interfaces. This is somewhat different from the result found in our previous simulations,13 where the O-C1-C2-C3 dihedral angle was 40% trans at the water/vapor interface and 60% trans at the CCl4/ water interface. However, those simulations were performed for a much lower surface coverage and the force field for the dihedral potential for the head group was different. These features can account for the observed difference. The probability of the gauche defect for next carbons oscillates, while for the last carbons it increases slightly. The average number of gauche defects in the chain for the monolayer at the CCl4/water and the water/vapor interface are 2.6 and 2.8, respectively. 3.3. Electrical Potential. The value of the interface (dipole) potential depends on the value of charges distributed in the system. In simulations the charge distribution depends on the force field used. However, even in the framework of the given force field the value of the potential is very sensitive to how the long-range electrostatic forces are treated in the simulation. In the present and previous simulations we use an Ewald summation technique to calculate the long-range Coulombic forces. However, since most of the interfacial systems are not periodic in one of the three dimensions (for instance the Z-dimension), conventional three-dimensional (3D) Ewald summation is not directly applicable. Instead, a more computationally expensive two-dimensional (2D) Ewald summation technique should be used.28-30 Another frequently used approach is to apply a conventional 3D Ewald summation method to a box that has a sufficiently large size in the z-direction, so that empty space between periodic images in the z-direction is present in the system.30-32 The empty space is created to avoid the artificial influence from these periodic images. For symmetric systems that have a zero dipole moment, using the 3D Ewald summation technique with the sufficiently large length of the simulation cell in the z-direction produces results that converge to the results obtained from simulations with a 2D Ewald sum. When asymmetric systems are simulated, i.e.,

Figure 6. Charge density profile for the SDS monolayer at the CCl4/ water interface with and without the dipole correction term (top panel). The electric field is depicted in the middle panel and the potential profile in the bottom. The dashed line is for the system with the correction term included, and the solid one is for the system without the correction in the Ewald summation.

systems with a net total dipole moment (in the z-direction), the 3D Ewald summation method has to be used with caution. As it was shown recently,21 to use the 3D Ewald method in this case, one has to add the surface term to the expression of electrostatic energy. As a result, we need to add extra contributions to the forces and these are given by the following formulas:21

Fx,i ) Fy,i ) 0 Fz,i ) -

qiMz  0V

(3)

where Mz is the z-component of the total dipole moment, qi is the charge of the i-molecule, and V is the volume. To find out how the presence of the extra force influences the structural and electrical properties of the system, we carried out simulations of SDS molecules at the CCl4/water interface with and without the correction term. The two simulations produced no substantial difference in the density of the liquids and surfactants. We expected that the major difference could appear in electric properties of the interface such as the surface potential. Thus, we calculated the electric field and the electric potential for our system which are displayed in Figure 6. The electric field is calculated using the formula

Ez(z) )

1 0

∫zz

1

dz′ Fq(z′)

(4)

and the potential is given by the formula

∆φ ) φ(z2) - φ(z1) ) -

∫zz

2

1

dz′ Ez(z′)

(5)

where Fq(z) is the charge density. The reference potential φ(z1) ) 0.0 V was chosen in the vacuum region (z < 0) far from the interface.

5306 J. Phys. Chem. B, Vol. 104, No. 22, 2000

Dominguez and Berkowitz how it is affected by the presence of the monolayer surfactant, we analyzed the structure of the interface in slightly more detail. For this purpose we followed the procedure used in the previous studies of liquid/liquid interfaces.33,34 According to this procedure, additional details on a surface width can be obtained from the measurement of a local surface width (wij). To obtain this quantity, we divided the XY-plane into N × N squares and considered at each instance all the molecules in the ijth box of cross section XY/N2. The largest N should be such that xS/N > ξb, where S is the area of the system and ξb is the bulk liquid correlation length (about 4-5 Å for most liquids far from the critical point). This way we do not consider the fluctuations that have a length scale of ξb. In each column we calculate the local width wij using the following expression:

wij(t) ) max[zH20(i,j)] - min[zCCl4(i,j)]

Figure 7. Snapshot of a configuration of SDS monolayer at the CCl4/ water interface. Black stick-balls represent SDS molecules, gray molecules on the top are CCl4, and the dark gray molecules at the bottom are water (bigger balls represent oxygen and smaller balls hydrogen).

where max[zH20(i,j)] is the position of the oxygen with the largest z-coordinate in the ij column and min[zCCl4(i,j)] is the position of the carbon (in CCl4) with the smallest z-coordinate in the ij column. Values of wij we get from simulations can be used to define the average width of the interface at a given time step and the average width over the trajectory, which are defined as

〈w(t)〉 )

〈w〉 )

1 N

Figure 8. Snapshot of a SDS monolayer at the vapor/water interface.

We observed that the electric fields are similar in both simulations. The difference becomes more pronounced for the surface potential. The potential difference across the water/SDS/ CCl4 interface when the correction is included is approximately 240 ( 35 mV, whereas when the correction term is omitted the difference is 160 ( 11 mV (see Figure 6). Since the electric field due to the correction term is relatively small (Mz/0V ≈ 0.02 V/Å) in our case, the deviation in the potential is not that large. For systems with a very small dipole moment, such as the one containing vapor/CCl4/water/vapor interfaces, the difference in the potential obtained from simulations with and without the correction term is negligible, as we observed in our simulations (not shown here). For systems such as water between two Pt(111) walls in the presence of a strong external electric field,21 the correction term can be large, leading to a large difference in the electric field and potential, depending on the inclusion of the extra term. 3.4. Surface Structure of the CCl4/Water Interface. We present snapshots from our simulations of the SDS monolayer at the CCl4/water and the water/vapor interfaces in Figures 7 and 8, respectively. As the figures show, the head groups of SDS molecules are well solvated by water, whereas when the organic liquid is present a considerable number of CCl4 molecules solvate the tails of the surfactant molecules. In order to study the width of the CCl4/water interface and to understand

(6)

1 N2

wij(t) ∑ i,j

1

∫wij(t) dt ∑ 2 i,j T

(7)

(8)

Results for the distributions of average widths 〈w(t)〉 are shown in Figure 9 for N ) 1, 2, 4, and 6 . We note that the distributions become more narrow and that the peaks of these distributions are shifted toward negative values of the z-coordinate as N increases. In particular the distribution for the system water/ SDS/CCl4 is shifted more (as N increases) compared to the pure water/CCl4 system. This indicates that water and CCl4 may not be in contact everywhere (i.e., there are regions where water and CCl4 are well separated; see Figure 7). Moreover, the fact that the width distribution for both systems seems to depend on the number of squares in the surface indicates that the interface is not flat, even in the absence of the monolayer. The distribution functions represented in Figure 9 also allow us to quantify the meaning of the surface roughness. Since the width of the interface depends on the resolution of our observation (number of squares), we will measure the width when the number of squares is the largest, but the fluctuations on a scale of bulk correlation are smoothed out; i.e., we look at the distribution for N ) 6 . From the distribution of the widths we obtained that at this resolution the water/CCl4 has a width of ∼3.3 Å, while in the presence of the surfactant the width is ∼8.4 Å. 3.5. Surface Tension. To find the values of the surface tensions, we calculated the quantity Γ, given by the following equation:

Γ ) Lz(〈Pn〉 - 〈Pt〉)

(9)

where 〈Pn〉 is the normal pressure, 〈Pt〉 is the average tangential pressure, and Lz is the length of the box in the z-direction. Since the calculation of the pressure involves the calculation of the virial, it is clear that the value of the pressure will depend on the inclusion of the surface term given by eq 3. Below we

SDS at Liquid/Liquid and Liquid/Vapor Interfaces

J. Phys. Chem. B, Vol. 104, No. 22, 2000 5307 of 72 and 26.9 mN/m.35,36 For the system vapor/water/SDS/ CCl4/vapor we obtained Γ ) 90.6 ( 5 mN/m when no dipole correction term was included. For the same system Γ became 100.7 ( 6 mN/m when the dipole correction term was used. For the simulation of the system vapor/water/CCl4/vapor the value of Γ was 131.8 ( 4 mN/m. Therefore from eq 11 we obtained that γw/o ) 41.6 ( 11 mN/m and γ˜ w/o ) 10.5 ( 13 mN/m. As we can see, the presence of surfactant at the interface reduced the value of the surface tension. It is interesting that the value of this reduction depends on the presence of the correction term. If we perform the same calculation using the value for Γ obtained from the calculation without the dipole correction term, we obtain that γ˜ w/o ) 0.4 ( 12 mN/m. Although the values of γ˜ w/o when calculated with and without the correction term are close to each other within the limits of the error bars, the value calculated when the correction dipole term was included makes more sense. 4. Conclusions We performed molecular dynamics computer simulations to study the structural properties of SDS monolayers at water/CCl4 and water/vapor interfaces. The area per surfactant molecule in our simulations was 45 Å2. We found that at this surface coverage tails of surfactant molecules are disordered and inclined with respect to the bilayer normal. The inclination is larger for SDS molecules at the water/vapor interface compared to the water/CCl4 interface and so is the disorder, although the difference is very small. We also measured the dipole potential and the surface tension of the interface. We showed that for these properties the inclusion of the dipole correction term into the calculation of long-ranged electrostatic forces is important if one wants to obtain good quantitative results.

Figure 9. Probability distribution for the width PN(〈w(t)〉) of the CCl4/ water interface. The distributions were calculated using different surface areas (N × N). Dashed lines are for the CCl4/water interface, and solid lines are for CCl4/water interface with the SDS monolayer.

Acknowledgment. Simulations were performed on a SGI Origin 2000 at the University of North Carolina and on a CrayT3D at the North Carolina Supercomputer Center. The work was supported by a grant from the Office of Naval Research. We are grateful to Dr. A. M. Smondyrev for his help and advice in performing the simulations and careful reading of the manuscript. References and Notes

describe how the inclusion of this term changes the result for the surface tension of the interface. When only one interface is present in the system, Γ is equal to the surface tension γ of this interface. When the simulation box contains several well-separated interfaces

Γ)

∑γi

(10)

In case of the vapor/water/CCl4/vapor system

Γ ) γv/w + γw/o + γo/v

(11)

where γv/w is the surface tension for the vapor/water interface, γw/o is the surface tension for the water/organic liquid interface, and γo/v is the surface tension for the organic liquid/vapor interface. In the presence of SDS molecules at the water/CCl4 interface the same formula can be employed; only the surface tension of the interface water/organic liquid γ˜ w/o will be changed (it will become smaller). To calculate this decrease, we need to find the values of the surface tensions for the vapor/water and CCl4/vapor interfaces. Values obtained from the simulations of slabs of these liquids are γv/w ) 68.8 ( 4 mN/m and γo/v ) 21.4 ( 3 mN/m, in good agreement with the experimental values

(1) Piasecki, D. A.; Wirth, M. J. J. Phys. Chem. 1993, 97, 7700. (2) Lu, J. R.; Marrocco, A.; Su, T. J.; Thomas, R. K.; Penfold, J. J. Colloid Interface Sci. 1993, 158, 303. (3) Tian, Y.; Umemura, J.; Takenaka, T. Langmuir 1988, 4, 1064. (4) Higgins, D. A.; Naujok, R. R.; Corn, R. M. Chem. Phys. Lett. 1993, 213, 485. (5) Conboy, J. C.; Messmer, M. C.; Richmond, G. J. Phys. Chem. B 1997, 101, 6724. (6) Conboy, J. C.; Messmer, M. C.; Richmond, G. J. Phys. Chem. 1996, 100, 7617. (7) Conboy, J. C.; Messmer, M. C.; Richmond, G. Langmuir 1998, 14, 6722. (8) Zhang, Z. H.; Tsuyumoto, I.; Kitamori, T.; Sawada, T. J. Phys. Chem. B 1998, 102, 10284. (9) Lyttle, D. J.; Lu, J. R.; Su, T. J.; Thomas, R. K.; Penfold, J. Langmuir 1995, 11, 1001. (10) Lu, J. R.; Hromadova, M.; Simister, E. A.; Thomas, R. K.; Penfold, J. J. Phys. Chem. 1994, 98, 11519. (11) Kuhn, H.; Rehage, H. J. Phys. Chem. B 1999, 103, 8493. (12) Michael, D.; Benjamin, I. J. Phys. Chem. B 1998, 102, 5145. (13) Schweighofer, K. J.; Essmann, U.; Berkowitz, M. J. Phys. Chem. B 1997, 101, 3793. (14) Schweighofer, K. J.; Essmann, U.; Berkowitz, M. J. Phys. Chem. B 1997, 101, 10775. (15) Tarek, M.; Tobias, D. J.; Klein, M. L. J. Phys. Chem. 1995, 99, 1393. (16) Pohorille, A.; Benjamin, I. J. Phys. Chem. 1993, 97, 2664. (17) Pohorille, A.; Benjamin, I. J. Chem. Phys. 1991, 94, 5599.

5308 J. Phys. Chem. B, Vol. 104, No. 22, 2000 (18) Gao, J.; Jorgensen, W. L. J. Phys. Chem. 1988, 92, 5813. (19) Dominguez, H.; Smondyrev, A. M.; Berkowitz, M. L. J. Phys. Chem. B 1999, 103, 9588. (20) Okamura, E.; Fukushima, N.; Hayashi, S. Langmuir 1999, 15, 3589. (21) In-Chul, Y.; Berkowitz, M. L. J. Chem. Phys. 1999, 111, 3155. (22) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular forces; Pullman, B., Ed.; Reidel: Dordrecht, 1981; p 311. (23) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagonax, G.; Profeta, S.; Weiner, P. J. Am. Chem. Soc. 1984, 106, 765. (24) Forester, T. R.; Smith, W. DL-POLY Package of Molecular Simulations; CCLRC, Daresbury Laboratory: Daresbury, Warrington, England, 1996. (25) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (26) Smith, E. R. Proc. R. Soc. London, A 1981, 375, 475.

Dominguez and Berkowitz (27) Egberts, E.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 3718. (28) Parry, D. E. Surf. Sci. 1975, 49, 433. (29) Heyes, D. M.; Barber, M.; Clarke, J. H. R. J. Chem. Soc. Faraday Trans. 2 1977, 73, 1485. (30) Spohr, E. J. Chem. Phys. 1997, 107, 6342. (31) Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385. (32) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574. (33) Benjamin, I. J. Chem. Phys. 1992, 97, 1432. (34) Michael, D.; Benjamin, I. J. Phys. Chem. 1995, 99, 1530. (35) Grubb, S. G.; Kim, M. W.; Rasing, Th.; Shen, Y. R. Langmuir 1988, 4, 452. (36) Tsujii, K. In Surface ActiVity, Principles, Phenomena and Applications; Academic Press: New York, 1998; p 4.