11720
J. Phys. Chem. 1996, 100, 11720-11725
Molecular Dynamics Simulations of the Liquid/Vapor Interface of SPC/E Water Ramona S. Taylor, Liem X. Dang, and Bruce C. Garrett* EnVironmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, Washington 99352 ReceiVed: February 29, 1996X
Molecular dynamics computer simulations have been used to explore the structural and dynamical properties of water’s liquid/vapor interface using the simple extended point charge (SPC/E) model. Comparisons to the existing experimental and simulation data suggest that the SPC/E potential energy function provides a semiquantitative description of this interface. The orientation of H2O molecules at the interface is found to be bimodal in nature. The self-diffusion constant of water is calculated to be larger at the surface than in the bulk.
I. Introduction The structure and dynamics of the liquid/vapor interface of water have been the subject of considerable experimental and computational interest in recent years.1-23 Much of this interest is due to the role that this interface plays in not only chemistry but also biology and atmospheric science. Measurements of average quantities such as the surface tension and surface potential have long been possible. However, only with recent advancements in the experimental techniques of second harmonic generation1,2,24 and sum frequency generation3,4 has it been possible to obtain detailed information about the orientation and environment of the individual interfacial H2O molecules. Over the past two decades, these properties have also been examined using both theory and simulations.5-23 Yet very little has been done in the way of examining the dynamical properties of this important interface. Comparisons of the dynamical properties of bulk and interfacial water have only been made for water at 298 K.7,22,23 Here, we examine these properties over a 100 K temperature range. In at least one proposed mechanism for the uptake of gaseous molecules by H2O droplets in the troposphere, the dynamics of the droplet surface is believed to play a role in the overall process.25 The incoming gaseous particle is incorporated into the droplet’s bulk via a mechanism in which the surface H2O molecules surround and subsequently transport the particle inwards. In those instances in which the mass accommodation of the incoming molecule involves a chemical reaction between the incoming gaseous particle and a H2O molecule, the replenishment of the surface with fresh H2O molecules may also be of extreme importance. Molecular dynamics (MD) computer simulations can be used to gain insight into these processes. The accuracy of any MD simulation, however, is limited by the accuracy of the interaction potential functions being used. Potential functions are available which reproduce many of the properties of bulk H2O.26-29 Yet, to a large degree, the ability of these potentials to reproduce the properties of the liquid/vapor interface of H2O remains undetermined. Although the extended simple point charge (SPC/E) potential26 does remarkably well at reproducing many of the bulk properties of H2O, at the liquid/vapor interface the environment of the individual H2O molecules is no longer bulklike. Each interfacial H2O molecule has fewer neighboring H2O molecules than it would if it were located in the bulk of the liquid. For example, the values of the point charges in the SPC/E model approximate the average dipole moment of a X
Abstract published in AdVance ACS Abstracts, June 15, 1996.
S0022-3654(96)00615-6 CCC: $12.00
single water molecule in bulk water, but it is uncertain if this charge distribution is accurate in the interface region. In this paper, we use the SPC/E potential26 in conjunction with molecular dynamics computer simulations to explore both the structural and dynamical properties of the liquid/vapor interface of H2O. While this paper was in preparation, a similar set of calculations by Alejandre and co-workers was published.15 This work also calculates the surface tension of water using the SPC/E interaction potential. They, however, have used the Ewald summation method to calculate the long-range interactions. Thus, our paper enables a comparison of how a property such as surface tension is affected by the use of Ewald sums. In addition, unlike Alejandre and co-workers,15 we have also calculated the orientation and dynamical properties of SPC/E water. Comparisons of the results of these simulations to existing experimental data suggest that the SPC/E potential26 adequately describes this interface. The orientation of H2O molecules varies as a function of distance from the interface. At the vapor side of the interface, the dipole moment of the H2O molecules is directed away from the bulk liquid, while at the liquid side it is directed slightly toward the bulk liquid. This structure persists over the entire temperature range studied. The self-diffusion constant of water is calculated to be larger at the surface than in the bulk. Also, the residence time for an individual H2O molecule at the surface is calculated to be on the order of 2 ps. The outline of this paper is as follows. The details of the molecular dynamics simulations are given in section II. The results of the density distribution analysis, surface tension calculations, orientational distributions, and dynamical analysis are discussed in section III. The conclusions can be found in section IV. II. Method The water potential employed in these calculations is the SPC/E water potential developed by Berendsen and co-workers.26 For each water molecule, the OH bond length is defined to be 1.0 Å and the HOH bond angle is defined to be 109.47°. This potential treats the intermolecular interactions via a combination of a Lennard-Jones 6-12 potential function and a three-site electrostatic potential:
[( ) ( ) ]
U(rij) ) 4OO
σOO rOO
12
-
σOO rOO
© 1996 American Chemical Society
6
3
3
+ ∑∑ i)1 j)1
qiqj rij
(1)
Liquid/Vapor Interface of SPC/E Water
J. Phys. Chem., Vol. 100, No. 28, 1996 11721
Figure 1. The simulation cell. A rectangular section of water is sandwiched between two sections of “vapor”. The z axis is perpendicular to the two free water surfaces.
where qi and qj are the charges centered on the individual atoms of different H2O molecules and the Lennard-Jones interaction is calculated only between oxygen atoms of neighboring H2O molecules. The Lennard-Jones and electrostatic parameters are σOO ) 3.166 Å, OO ) 0.1554 kcal‚mol-1, qH ) +0.4238e, and qO ) -2.0qH. All of the molecular dynamics simulations were run using a modified version of the AMBER 4.0 suite of programs.30 AMBER 4.0 makes use of a Verlet integrator31 to solve the equations of motion. A step size of 0.001 ps was employed. Heating was accomplished via the Berendsen scheme32 with a coupling constant of 0.2. The OH and HH bond lengths and the intermediate HOH bond angle were constrained via the SHAKE algorithm.33 All simulations utilized a cutoff distance of 9 Å. For each temperature, the initial water configuration consisted of a rectangular section of H2O sandwiched between two sections of “vapor” as shown in Figure 1. Periodic boundary conditions are applied in all three directions. The initial lengths of the rectangle are determined by equilibrating a bulk box of H2O for a period of 200 ps in a NPT ensemble with a pressure of 1 atm.34 This results in different box lengths for each of the temperatures studied and for each box size studied. Once equilibrated, 50 Å of vapor is added to the box length in the z direction (as shown in Figure 1) and this new system is equilibrated for another 100 ps in a NVT ensemble with no pressure monitoring. Thus, the z direction becomes the direction perpendicular to the two newly created H2O surfaces. Simulations were run in this manner at 268, 283, 298, 323, 348, and 373 K. These simulations each consisted of 526 water molecules. To test the effects of the simulation size, simulations with 1052 water molecules were also run at 298 K. The box lengths in the x and y directions are 24.94, 24.94, 24.99, 25.16, 25.25, 25.44, and 31.58 Å for the simulations done at 268, 283, 298 (526 H2O molecules), 323, 348, 373, and 298 K (1052 H2O molecules), respectively. At each temperature, the system is equilibrated for 100 ps before any analysis is performed. For the density profiles and the orientational information, data was collected and averaged over a period of 75 ps. For the surface tension measurements, data was collected in over at least five 50 ps simulations. For the diffusional analysis, data was collected in two 50 ps simulations. III. Results and Discussion A side view of the H2O liquid/vapor interface equilibrated at 298 K is shown in Figure 2a. Neither of the two interfaces is smooth. Each may instead be best described as a landscape containing many valleys and hills. The region over which this roughness persists defines the surface thickness. In this region, the density of the water changes from that of the bulk (≈1
Figure 2. A snapshot of H2O at 298 K. The lighter gray H2O molecules are surface molecules as determined from the density profile (see text). (a, top) Side view with the z axis perpendicular to the two free surfaces. (b, bottom) Top view of the same system shown in (a). Here the z axis is perpendicular to the plane of the page.
g/cm3) to that of the vapor (≈0 g/cm3). The surface H2O molecules (light gray) were chosen as such from the density profile distribution described below. A top view of this same system is given in Figure 2b. From this picture, one can see that this is not a rigorous definition of a surface molecule and that some of the surface is covered with bulk H2O molecules. This problem will be addressed again in detail below. The density profile of H2O at 298 K is shown in Figure 3. This profile is constructed by dividing the simulation cell shown in Figure 1 into 1 Å thick layers along the z axis and subsequently calculating the H2O density in each layer. The resulting density profile exhibits two interfaces with the plateau region having a density consistent with that of bulk H2O. The density profile can be fit with a hyperbolic tangent function of the following form:6,11
[ ]
(z-z0) 1 1 F(z) ) (FL + FV) - (FL - FV) tanh 2 2 d
(2)
11722 J. Phys. Chem., Vol. 100, No. 28, 1996
Taylor et al.
Figure 3. Density profile of 526 water molecules equilibrated at 298 K. The open circles correspond to the simulation data and the solid line is the fit to this data obtained with the hyperbolic tangent function given in eq 2.
TABLE 1: Results of Fitting the z-Dependent Density Profiles with the Hyperbolic Tangent Function Shown in Eq 2a T (K)
FL (g cm-3)
zL (Å)
zR (Å)
t (Å)
γ (dyn cm-1)
268 283 298 298b 323 348 373
1.015 1.009 1.003 1.002 0.985 0.962 0.940
25.55 25.50 25.28 25.37 25.42 25.27 25.56
50.43 50.54 50.41 56.96 50.66 50.90 51.43
2.67 2.84 3.32 3.57 3.42 4.21 4.50
72((3) 68((2) 65((2) 66((2) 58((2) 51((1) 47((1)
Figure 4. Temperature dependence of the “10-90” interface thickness, t. The solid diamonds are the experimental ellipsometry data from ref 33. The solid squares are the results of the MD simulations using the CC H2O potential (ref 11). The solid circles are the results of this work and the open circles are the results of similar simulations done with the SPC/E potential (ref 15). In all cases the lines are only meant to guide the eye.
a T is the temperature of the system; F is the density of the bulk L liquid phase, zL and zR are the positions of the two individual Gibb’s surfaces, and t is the interface thickness. Also shown are the results of the surface tension (γ) calculations with the associated statistical uncertainties shown in parentheses. Unless otherwise indicated, all simulations were performed with 526 H2O molecules. b These simulations performed with 1056 H2O molecules.
where FL and FV are the densities corresponding to the liquid and vapor phases, z0 is the position of the Gibb’s dividing surface (where F(z) ) 1/2(FL + FV)), and d is the thickness parameter for the interface. The “10-90” thickness (t) is related to d via the following equation: t ) 2.1972d. The “10-90” thickness is the thickness over which the density of water changes from 0.9FL to 0.1FL assuming FV is equal to zero. The parameters resulting from the fitting of the density profiles and the corresponding “10-90” thickness are given in Table 1. (The calculated surface tension values are discussed in the next section.) In all cases, FV was found to be equal to zero and thus is not included in the table. The interface thickness increases with increasing temperature. The results of the simulations with N ) 1052 at 298 K are also given in the table. These results are the same within the error of the calculations as those found for the simulations done at 298 K with 526 H2O molecules. The interface thickness, t, has experimentally been measured via both ellipsometry35 and X-ray reflectivity36,37 and is plotted along with the results of these calculations and the results of similar calculations done by Matsumoto and Kataoka11 using the Carravetta-Clementi (CC) H2O potential38 in Figure 4. The comparison between the experimental data, these, and the Matsumoto and Kataoka11 simulations is not favorable. Both potentials underestimate the interface thickness. This is attributable to the lack of capillary waves in the simulations.11,18 Capillary waves are inherent in the experimental setup and are theoretically predicted to result in an expansion of the interface region. a. Surface Tension. The surface tension, or surface free energy (γ), is a measure of the degree of structure possessed by the interface or the desire of the molecules to be solvated.
Figure 5. Surface tension of H2O. The closed circles are the experimental data. The open circles are this work. The error bars represent one standard deviation of the simulation data. The ×’s are the results of MD simulations performed with the SPC/E potential and the Ewald summation method (ref 15).
Using the theories of Kirkwood and Buff,39-41 the surface tension of SPC/E water can be determined in a MD simulation using the following equation:
γ)
{ 〈∑ ∑ ( )( ) ∂VRβ
1 1
2 2A
i,j R,β
1
∂rRβ rRβ
〉}
(b r ij‚b r Rβ - 3zijzRβ)
(3)
where A is the area of the liquid/vapor interface, VRβ is the potential energy calculated between sites R and β on molecules i and j, rRβ and zRβ are the bond distance and the distance in the z direction between sites R and β, and rij and zij are the corresponding quantities calculated between the centers of mass of molecules i and j. The surface tension values of SPC/E water (open circles) calculated as a function of temperature are given in Table 1 and shown in Figure 5 along with the corresponding experimental data (closed circles).42 The error bars are one standard deviation of the five 50 ps trajectories. Given the simplicity with which the SPC/E potential treats the H2O molecule, the agreement between the MD simulations and the experimental data is quite good. Also shown in Figure 5 is the surface tension data from a recent paper by Alejandre, Tildesley, and Chapela.15 These authors also performed molecular dynamics simulations under similar conditions as ours with the SPC/E H2O potential. They, however, have calculated the H2O-H2O Coulombic interactions via the Ewald summation technique. By including the long-
Liquid/Vapor Interface of SPC/E Water
J. Phys. Chem., Vol. 100, No. 28, 1996 11723
Figure 6. Definition of the variables used in determining the orientation of the H2O molecules at the liquid/vapor interface. θ is defined as the angle between the surface normal and the H2O dipole vector, and φ is defined as the angle between the OH bond vector and the surface normal.
range Coulomb interactions in this manner, they obtain excellent agreement with the experimental data. The differences between our and their results arise from our neglect of those H2O-H2O interactions which occur beyond the 9 Å cutoff distance we are employing. It is encouraging, however, that the results of the two sets of simulations agree as well as they do. In the future it will be interesting to see how other interfacial properties, such as the dynamical properties described below, are affected by the inclusion of these long-range Coulomb interactions. b. Orientational Distributions. The orientational distribution of the H2O molecules at the liquid/vapor plays a large role in determining the reactivity of the interface. Recently, experimental techniques such as second harmonic generation (SHG) and sum frequency generation (SFG) have been used to probe the molecular structure of the liquid/vapor interface of H2O. Using SHG, Eisenthal and co-workers found that the dipole moment of H2O is slightly tilted inwards toward the bulk liquid.1,2 Using SFG, Shen and co-workers subsequently determined that 25% of the surface H2O molecules have a free OH bond which is directed out of the surface.3,4 To examine this orientational ordering in our simulations, it is necessary to define two angles between an individual H2O molecule and the surface normal. As shown in Figure 6, the first of these angles, designated θ, is defined as the angle between the H2O dipole vector and the surface normal. This angle can vary from 0° to 180° with an angle of 90° corresponding to the dipole vector lying in the plane of the interface. The second angle, φ, is defined as the angle between the OH bond and the surface normal. For each H2O molecule, two of these angles exist. The simulation system is divided up into 1 Å thick slices as described above for the density profile calculations. Thus, the orientation of the H2O molecules can be determined as a function of position. The orientational distributions of the H2O dipole vector for H2O equilibrated at 268, 298, and 373 K are given in Figure 7. In the liquid/vapor interfacial region, two distinct orientations of H2O molecules exist. On the vapor side of the interface, the most probable orientation is that with the H2O dipole directed out of the liquid at an angle of 74° relative to the surface normal. On the liquid side of the interface, however, the H2O molecules are oriented such that their dipole is lying in the plane of the interface. As expected, no orientation preference is found for the H2O molecules in the bulk portion of the simulation cell. This z-dependent shift in the dipole orientation peak has also been seen in the MD simulations of Matsumoto and Kataoka.11 In contrast to their work, however, we find that this structure is independent of temperature. A similar orientational distribution is found for each of the simulations performed. As shown in Figure 7, however, the intensity of the peaks decreases at the higher temperatures. The individual H2Os are energetically hotter and thus are less likely to be caught in one predominate configuration. The orientation distribution of H2O molecules at the liquid side of the interface persists over a longer distance
Figure 7. Orientation distribution of the H2O dipole vector relative to the surface normal as a function of temperature. The dual nature of the interface is shown by the various curves. The solid curve corresponds to the orientation of H2O molecules in the bulk of the liquid. The longer dashed curve corresponds to the orientation of molecules 3 Å from the surface. The shorter dashed curve corresponds to molecules 2 Å from the surface, and the dotted curve corresponds to the orientation distribution of the H2O molecules at the surface.
and thus likely provides the dominant portion of the signal in experiments such as the SHG work of Eisenthal and coworkers.1,2 The availability of free OH bonds at the liquid/vapor interface of H2O will affect the interface’s reactivity toward adsorbing particles. As seen in Figure 8 for H2O equilibrated at 268, 298, and 373 K, the peaks corresponding to the orientational distribution of the OH bond vectors shift toward higher angles relative to the surface normal as a function of z. At the vapor side of the interface, the predominate structure is one in which an H is pointed out of the surface. However, on the liquid side of the interface the hydrogens lie within the plane of the interface or are directed into the bulk with an angle of 94° relative to the surface normal. This ordering is again found to persist over the entire temperature range studied. c. Dynamical Properties. As indicated above, accurately defining the surface in a MD simulation of a liquid is not straightforward. Up until now, we have been able to divide the simulation cell into layers and calculate average quantities for the various structural properties of interest. However, as seen in Figure 9 where the outermost H2O molecules are depicted in light gray, the surface is not a smooth entity and thus cannot be defined by either a single layer or a set of parallel layers. The peaks and valleys in the surface often span many angstroms, and thus many layers. To differentiate between the
11724 J. Phys. Chem., Vol. 100, No. 28, 1996
Taylor et al.
Figure 10. Self-diffusion coefficient of both the surface and bulk H2O molecules as a function of temperature. The solid diamonds correspond to the self-diffusion coefficient of the surface H2Os in the plane of the H2O interface. The solid circles represent the diffusion coefficients of the “bulk” H2O molecules. The open circles correspond to the experimental bulk H2O diffusion coefficients measured via the proton spin echo technique (ref 42).
Figure 8. Orientation distribution of H2O’s OH bond relative to the surface normal as a function of temperature. The definition of the various curves is the same as that given for Figure 7.
cell shown in Figure 1 into columns which are 3 Å by 3 Å in the x and y directions and run the entire length of the simulation cell in the z direction. The outermost H2O molecules in each column (this would be two H2Ossone at each surface) are designated as surface molecules and the remaining molecules are designated as bulk molecules. Even this definition is not rigid enough. A gray area remains as to whether molecules just below the surface are more bulklike or more surfacelike. However, it does provide a cleaner definition of the surface molecules than that used in the calculation of the structural properties. The self-diffusion coefficient of the H2O molecules at the surface and in the bulk can be calculated from the Einstein relation:43
D)
Figure 9. Snapshot of one surface of H2O equilibrated at 298 K. The light gray molecules represent the outermost H2Os. The penetration depth of the surface H2O is clearly visible.
dynamics of the surface molecules and the bulk molecules, a new way of defining or differentiating the surface H2O molecules is needed. This was done by dividing the simulation
〈|ri(t;t*) - ri(t0;t*)|2〉 2f∆t
(4)
where D is the diffusion coefficient, ri(t;t*) is the position of particle i at time t, ri(t0;t*) is the position of particle i at time 0, ∆t is the elapsed time between t and t0 assuming that the molecule has not left the interface for longer than t*, and f is the number of directional degrees of freedom. In the bulk, f ) 3 and the sum |ri(t;t*) - ri(t0;t*)| includes the x, y, and z contributions. At the interface f ) 2 and the sum |ri(t;t*) ri(t0;t*)| only includes the x and y contributions. D is calculated from a least-squares fit to the slope of the plot of 〈|ri(t;t*) ri(t0;t*)|2〉 vs ∆t. For a surface molecule, only the diffusion coefficient for movement within the surface plane can be calculated since if a molecule moves perpendicular to the surface plane in a downward direction it becomes a bulk molecule. If a molecule leaves the surface for less than 2 ps it is considered to remain a surface molecule. The choice of t* ) 2 ps will be discussed below. Each time the molecule changes from being a surface molecule to being a bulk molecule, ri(t0) and ∆t are reinitialized. This requires a longer total simulation time so that good statistics can be obtained for those molecules which stay at the surface for periods of time larger than 15 ps. The diffusion coefficients calculated in this way as a function of temperature are shown in Figure 10. At all temperatures, diffusion is faster at the surface than in the bulk. This is in agreement with the previous findings7,22,23 for water at 298 K. At the liquid/vapor interface, each H2O molecule has on average two less hydrogen bonds than it would have in the bulk. This breakdown of the hydrogen-bonding network allows the interfacial water molecules to move more freely than those in the
Liquid/Vapor Interface of SPC/E Water
J. Phys. Chem., Vol. 100, No. 28, 1996 11725
bulk, and thus their diffusion coefficient is larger. The diffusion coefficient of bulk water is slightly larger than that found experimentally.26,44 This difference, however, is due to the gray area described above and not to the SPC/E potential function. Diffusion coefficients calculated from purely bulk water using this potential agree quite well with the experimental data.26 Water molecules are constantly moving from the bulk to the interface and back. It is not clear how long after a molecule leaves the liquid/vapor interface that it loses its surface character. The probability that a H2O molecule will remain at the interface for a period of time greater than (t - t0) can be calculated from the following correlation function:45,46
n(t - t0) ) 〈∑Pi(t,t0;t*)〉
(5)
i
where Pi(t,t0;t*) ) 1 if H2O molecule i is a interfacial molecule at both times t and t0, and if the water molecule has left the interface for no longer than t* in the interim. Otherwise, Pi(t,t0;t*) ) 0. The average lifetime, or residence time, τR, of a H2O molecule at the interface is defined as the exponential decay time of the above correlation function (eq 5). τR is also dependent upon t*, the amount of time that the molecule leaves the interface and is still thought to be a surface molecule. However, for a molecule at the H2O surface at 298 K, τR is insensitive to the value of t*. Within the error of the calculation, varying the value of t* from 0 to 2 ps does not affect the resulting residence time of 2 ps. Thus, if a H2O molecule leaves the interface for longer than this average lifetime, τR, the probability of it retaining its interfacial character is low and it consequently is labeled as a bulk H2O molecule. IV. Conclusions Molecular dynamics computer simulations have been used to explore the structural and dynamical properties of the liquid/ vapor interface of SPC/E water. The calculated surface tension of H2O compares favorably with the corresponding experimental measurements. The orientation of the liquid H2O molecules changes as a function of distance from the vapor. The dipoles of the H2O molecules at the vapor side of the interface are directed out of the surface at an average angle of 74° relative to the surface normal and the dipole of the molecules at the liquid side of the interface lie parallel to the surface plane. At the vapor side of the interface, one water hydrogen points on average straight out of the surface, while both are pointed slightly toward the bulk on the liquid side of the interface. The self-diffusion coefficient of the interfacial H2O molecules is larger than that of the bulk molecules. At 298 K, the residence time for a H2O molecule at the interface is on the order of 2 ps. The mechanisms of Worsnop and Davidovits25 rely on the ability of the “interfacial” water molecules to encompass the incoming molecule and draw it into the bulk of the liquid as opposed to a cavity forming under the molecule into which the molecule can fall into the bulk liquid. Understanding the dynamics of the surface molecules may be of help in refining these mechanisms. The effect of adding an adsorbate to the interface may also be of interest. Simulations using the SPC/E water potential which are aimed at better understanding the uptake of gaseous particles by H2O droplets are underway. Acknowledgment. We thank Dr. Doug Ray and Dr. Greg Schenter for insightful discussions. This work was performed under the auspices of the Division of Chemical Science, U.S.
Department of Energy, under contract DE-AC06-76RLO 1830 with Battelle Memorial Institute, which operates the Pacific Northwest National Laboratory. R.S.T. also acknowledges a postdoctoral fellowship from Associated Western Universities, Inc. under grant No. DE-ACO6-76RLO 1830 with the U.S. Department of Energy. References and Notes (1) Goh, M. C.; Hicks, J. M.; Kemnitz, K.; Pinto, G. R.; Bhattacharyya, K.; Eisenthal, K. B.; Heinz, T. F. J. Phys. Chem. 1988, 92, 5074. (2) Goh, M. C.; Eisenthal, K. B. Chem. Phys. Lett. 1989, 157, 101. (3) Du, Q.; Superfine, R.; Freysz, E.; Shen, Y. R. Phys. ReV. Lett. 1993, 70, 2313. (4) Du, Q.; Freysz, E.; Shen, Y. R. Science 1994, 264, 826. (5) Stillenger, F. H., Jr.; Ben-Naim, A. J. Chem. Phys. 1967, 47, 4431. (6) Townsend, R. M.; Gryko, J.; Rice, S. A. J. Chem. Phys. 1985, 82, 4391. (7) Townsend, R. M.; Rice, S. A. J. Chem. Phys. 1991, 94, 2207. (8) Wilson, M. A.; Pohorille, A.; Pratt, L. R. J. Phys. Chem. 1987, 91, 4873. (9) Wilson, M. A.; Pohorille, A.; Pratt, L. R. J. Chem. Phys. 1988, 88, 3281. (10) Wilson, M. A.; Pohorille, A.; Pratt, L. R. J. Chem. Phys. 1989, 90, 5211. (11) Matsumoto, M.; Kataoka, Y. J. Chem. Phys. 1988, 88, 3233. (12) Motakabbir, K. A.; Berkowitz, M. L. Chem. Phys. Lett. 1991, 176, 61. (13) Duncan-Hewitt, W. C. Langmuir 1991, 7, 1229. (14) Lie, G. C.; Grigoras, S.; Dang, L. X.; Yang, D.-Y.; McLean, A. D. J. Chem. Phys. 1993, 99, 3933. (15) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. J. Chem. Phys. 1995, 102, 4574. (16) Benjamin, I. Phys. ReV. Lett. 1994, 73, 2083. (17) Barraclough, C. G.; McTigue, P. T.; Ng, Y. L. J. Electroanal. Chem. 1992, 329, 9. (18) Weeks, J. D. J. Chem. Phys. 1977, 67, 3106. (19) Croxton, C. A. Phys. Lett. 1979, 74A, 325. (20) Croxton, C. A. Physica 1981, 106A, 239. (21) Besseling, N. A. M.; Lyklema, J. J. Phys. Chem. 1994, 98, 11610. (22) Zhu, S.-B.; Fillingim, T. G.; Robinson, G. W. J. Phys. Chem. 1991, 95, 1002. (23) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1988, 89, 5044. (24) Eisenthal, K. B. Annu. ReV. Phys. Chem. 1992, 43, 627. (25) Davidovits, P.; Jayne, J. T.; Duan, S. X.; Worsnop, D. R.; Zahniser, M. S.; Kolb, C. E. J. Phys. Chem. 1991, 95, 6337. (26) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (27) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (28) Dang, L. X. J. Chem. Phys. 1992, 97, 2659. (29) Watanabe, K.; Klein, M. L. Chem. Phys. 1989, 131, 157. (30) Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Seibel, G. L.; Singh, U. C.; Weiner, P.; Kollman, P. AMBER 4.0, University of California, San Francisco, 1991. (31) Verlet, L. Phys. ReV. 1967, 159, 98. (32) Berendsen, H. J. C.; Postama, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (33) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (34) Andersen, H. C. J. Chem. Phys. 1980, 72, 2384. (35) Kinosita, K.; Yokota, H. J. Phys. Soc. Jpn. 1965, 20, 1086. (36) Braslau, A.; Deutsch, M.; Pershan, P. S.; Weiss, A. H.; Als-Nielsen, J.; Bohr, J. Phys. ReV. Lett. 1985, 54, 114. (37) Braslau, A.; Pershan, P. S.; Swislow, G.; Ocko, B. M.; Als-Nielsen, J. Phys. ReV. 1988, A38, 2457. (38) Carravetta, V.; Clementi, E. J. Chem. Phys. 1984, 81, 2646. (39) Fowler, R. H. Proc. R. Soc. 1937, A159, 229. (40) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1949, 17, 338. (41) Salomons, E.; Mareschal, M. J. Phys.: Condens. Matter 1991, 3, 3645. (42) MacRitchie, F. Chemistry at Interfaces; Academic Press, Inc.: San Diego, CA, 1990; p 269. (43) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1987; p 60. (44) Krynicki, K.; Green, C. D.; Sawyer, D. W. Faraday Discuss. Chem. Soc. 1979, 66, 199. (45) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071. (46) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757.
JP960615B
