J . Phys. Chem. 1987, 91, 4873-4878
0.
Intermolecular potential
c
'-1
I
?
i
1;
-2000.
A
f I
I
/
1000.
0.
Temperature (Kelvin) Figure 10. Intermolecular potential energy of the water dimer as a function of temperature. Symbols are as follows: (0)quantum simulation, with error bars; (A)classical simulation; (X) quantum result from Wallqvist and Berne;'O (+) classical result from same authors. The 0 K value is taken from a random walk simulation of the time-dependent
4873
degrees of freed~m,*~.~O they are likely to be inadequate for intramolecular vibrations. Figures 9 and 10 give the intramolecular and intermolecular contributions to the potential energy of the water dimer. There are large differences between the quantum and classical results for the intramolecular potential energy, the former being almost independent of temperature and close to the sum of the independent oscillator ground-state energies. Deviations from the ground-state energy that are apparent as the temperature is increased are largely due to excitation of the bending modes. So far as the intermolecular potential is concerned, the quantum and classical results converge above 300 K. Also shown in Figure 10 are results of path integral Monte Carlo calculations on the water dimer.30 The intermolecular potential function used by Wallqvist and Berne was the central force model of Stillinger and Rahman, although the intramolecular potential was that used here.27 Using the central force model in classical simulations gives an intermolecular potential energy which is about 4% higher than our classical value at 100 K, and similar differences are expected for the quantum result. The differences between quantum and classical simulation results, based on both quantum simulations, are about 10% at 100 K and 5% at 200 K.
structure is evident in the functions goH and gHH, corresponding to the increased stabilization of the hydrogen-bonded configuration. Only at temperatures below 300 K do differences between the classical and quantum calculations become apparent. From the intermolecular and intramolecular distribution functions we can conclude that, whereas classical simulations of liquid water35should give good results for properties associated with intermolecular
Summary We have demonstrated that the quantum Monte Carlo method, based on an analogy between the time-dependent Schrodinger equation and the diffusion equation with source and sink terms, can be used to solve the Bloch equation for the equilibrium density matrix. In this paper we have developed the method by considering various simple problems, comparing the calculations with both analytic and other simulation results. The method appears to be capable of giving accurate information on problems in quantum statistical mechanics.
(35) Watts, R. 0.;McGee, I. J. Liquid State Chemicd Physics; WileyInterscience: New York, 1976.
Acknowledgment. The work described in this paper was supported by a grant of computer time, on a Cyber 205, by the Institute of Physical Sciences, CSIRO, Canberra.
Schrodinger equation.*'
Molecular Dynamics of the Water Liquid-Vapor Interface Michael A. Wilson, Andrew Pohorille, Department of Chemistry, University of California, Berkeley, California 94720
and Lawrence R. Pratt* Chemical and Laser Sciences Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 (Received: October 22, 1986)
The results of molecular dynamics calculations on the equilibrium interface between liquid water and its vapor at 325 K are presented. For the TIP4P model of water intermolecular pair potentials, the average surface dipole density points from the vapor to the liquid. The most common orientations of water molecules have the C, molecular axis roughly parallel to the interface. The distributions are quite broad and therefore compatible with the intermolecular correlations characteristic of bulk liquid water. All near-neighbor pairs in the outermost interfacial layers are hydrogen bonded according to the common definition adopted here. The orientational preferences of water molecules near a free surface differ from those near rigidly planar walls which can be interpreted in terms of patterns found in hexagonal ice 1. The mean electric field in the interfacial region is parallel to the mean polarization which indicates that attention cannot be limited to dipolar charge distributions in macroscopic descriptions of the electrical properties of this interface. The value of the surface tension obtained is 132 f 46 dyn/cm, significantly different from the value for experimental water of 68 dyn/cm at 325 K.
Introduction This report presents the results of molecular dynamics calculations on an interface between liquid water and its vapor. Results of particular interest include the nature of ordering of water 0022-3654/87/2091-4873$01.50/0
molecules in an interfacial region and the adequacy of intermolecular Potential energy n-"M for treatment of surfaces. The surfaces of chemically complex aqueous systems are likely to be intensely studied by computer simulation methods over the 0 1987 American Chemical Society
4814
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987
next several years. This is due to the chemical importance of aqueous interfaces and to rapid improvement of computational capabilities which make well-controlled studies of these complicated systems feasible. A first step in such a direction is study of the interface between pure water liquid and vapor. Currently, computer experimental data on the free surfaces of pure liquid water are sparse. Thus, an objective of this work is to provide for this simple system information which might be used to better understand results of computer experiments on more complicated systems. The interfacial ordering of water molecules has been the subject of several previous theoretical contributions. Stillinger and Ben-Naiml and Croxton* each discussed the orientational preferences of water molecules near a liquid-vapor surface and arrived at predictions which differ with respect to the direction of the surface polarization. This issue has not been resolved by laboratory experiments. Preliminary computer simulation results on the water liquid-vapor interface have been presented b e f ~ r e . ~ Those .~ calculations, however, do not establish a clear view of the orientational ordering of water molecules near a fluid interface, nor do they provide benchmark data for comparison to results of calculations on more complicated aqueous surfaces. Another type of water surface has been studied by analyzing the minimum energy configurations of clusters of water molecule^.^ Considerably more computer experimental data are available for liquid water abutting a rigidly planar ~ a l l . ~ 3 Comparison ~-~ of the present results with results of previous calculations8 provides information about the differences in the molecular structures of free and constrained aqueous surfaces. Additionally, comparison of the present results with measured properties of water provides information about the fidelity of available intermolecular potential energy models to aqueous surfaces. In particular, we estimate the surface tension of TIP4P9 water at 325 K and compare the computed value with the surface tension of experimental water at the temperature. The next section documents the essential aspects of the molecular dynamic calculations. Subsequent sections present the results and a concluding discussion.
Methods The initial configuration of 342 water molecules in a lamellar arrangement was constructed from the final configuration of another molecular dynamics calculation on an aqueous interface.I0 The dimensions of the simulation cell were 21.71 X 21.71 A X 65.1 A. The TIP4P model potential was used, and the intermolecular interactions were smoothly truncated by a cubic spline modification between 6.5 and 7.0 A as discussed by Andersen." Periodic boundary conditions were adopted in all three spatial dimensions. The equations of motion were solved with the Verlet algorithmI2 and SHAKEI3 using a time step of 0.005 ps. A trajectory of 40 ps was constructed during which intermittent ad( 1 ) Stillinger, F. H.; Ben-Naim, A. J . Chem. Phys. 1967, 47, 4431. (2) Croxton, C. A. Statistical Mechanics of the Liquid Surface; Wiley: New York, 1980; Chapter 7; Phys. Lett. 1979, 74A, 325; Physica A (Amsterdam) 1981, lObA, 239. ( 3 ) Townsend, R. M.; Gryko, J.; Rice, S. A. J. Chem. Phys. 1985, 82, 4391. (4) Lee, C. Y.; Scott, H. L. J . Chem. Phys. 1980, 73, 4591. ( 5 ) Weber, T. A,; Stillinger, F. H. J . Phys. Chem. 1938, 87, 4277. (6) Jonsson, B. Chem. Phys. Lett. 1980, 82, 520. (7) Marchesi, M. Chem. Phys. Lett. 1983, 97, 224. (8) Lee, C. Y. McCammon, J. A.; Rossky, P. J. J . Chem. Phys. 1984,80, 4448. (9) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., J. Chem. Phys. 1983, 79, 926. (10) Wilson, M. A.; Pohorille, A.; Pratt, L. R., unpublished molecular dynamics results on Na' ion motion near a water liquid-vapor interface. For related calculations on a bulk solution see: Wilson, M. A. Pohorille, A.; Pratt, L. R. J. Chem. Phys. 1985.83, 5832. ( 1 I ) Andrea, T. A.; Swope, W. C.; Andersen, H . C. J. Chem. Phys. 1983, 79, 4571. (12) Verlet, L. Phys. Reu. 1967, 159, 98. (13) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23. 321
Wilson et al. justment of atomic velocities brought the system to a temperature of about 320 K. The system was then aged for an additional 30 ps. Finally, a trajectory of 140 ps was computed. The average temperature was 325 K. The present analysis is based on the 140 configurations collected by storing the configuration at the end of each 1 .O-ps segment of the trajectory. Statistical uncertainties were estimated on the basis of the assumption that each configuration represents a statistically independent observation. A technical difficulty encountered in calculations such as this one is that of aging the system long enough that the correct stationary state of phase equilibrium is achieved. For this system a close estimate of the densities of the coexisting phases could be made prior to starting the calculation and then used in initiating the calculation. The experimental volume of the liquid water in equilibrium with its vapor is 0.99 g/cm3 at 325 K.I4 Since the TIP4P potential is adjusted to fit the experimental pressure of bulk liquid at thermodynamic states similar to this one, we expect that it will predict reasonable densities for the homogeneous liquid portions of the simulated system under discussion. This liquid density is only slightly different from those most commonly studied in simulations of bulk aqueous phases. The density of saturated g/cm3 at the temperature of the calc~lation.'~ steam is 9.05 X If the value 21.71 X 21.71 X (65.1 - 21.71) A3 is taken as an estimate of the volume of the system which can be ascribed to the vapor, then we expect an average number of molecules in the vapor of 0.06. For almost all of the calculation, essentially none of the 342 molecules would be expected to be identifiable as "vapor". No significant transport of material across the interface is expected. Another concern about calculations such as this one is the necessity of treating sufficiently large systems that a macroscopic limit is practically achieved for all properties of interest. As was pointed out above, the density of TIP4P water under the conditions of the present calculation is almost identical with that of experimental water. Thus, the lateral dimension of the system corresponds closely to seven molecular spacings in bulk liquid water near its triple point. For similarly sized lamellar systems of noble gases at liquid-vapor coexistence near their triple point, it has been estimated that system size corrections should lower the computed surface tension by approximately lO%.I5 This is smaller than our statistical uncertainty in estimating the surface tension. Therefore, it has been assumed that the system is a sufficiently large system, and calculations on other sized systems have not been undertaken. Subsequent work which explicitly studies the dependence of surface properties of chemical interest on system size would be helpful. A comment on the decision to truncate the intermolecular interactions rather than employ an Ewald treatment of long-ranged forces is warranted in closing this section on methods. On the basis of currently available computer experimental data on bulk systems it is likely that the most satisfactory treatment of Coulomb range interactions involves an Ewald summation.l6-I9 The differences between results obtained with prudently truncated interactions and those obtained by inclusion of Coulomb range interactions are subtle, and, generally, only of major importance for explicit study of dielectric properties."-I9 An Ewald treatment of longranged forces is considerably more time consuming. Thus, for studies like this one which are principally limited by statistical precision and do not inquire directly into dielectric properties, truncation of intermolecular interactions while treating sufficiently (14) Keenan, J . H.; Keyes, F. G.; Hill, P. G.: Moore, 3. G. Steam Tables; Wiley: New York, 1969. (15) Weeks, J. D. J. Chem. Phys. 1977, 67, 3106. (16) Or related developments which take explicit account of dielectric behavior of the region exterior to the computer simulation sample, such as ref 19 and preceding work cited there. (17) The Problem of Long-Range Forces in the Computer Simulation of Condensed Media; NRCC Proc. No. 9; Ceperley, D., Ed.;(National Technical Inoformation Service: Springfield, VA, 1980). (18) De Leeuw, S. W.; Perram, J . W.; Smith, E. R. Annu. Reu. Phys. Chem. 1986, 37, 245. (19) Neumann, M. J. Chem. Phys. 1985,82,5663; J . Chem. Phys. 1986, 85, 1561.
Molecular Dynamics of the Water Liquid-Vapor Interface le2
-
The Journal of Physical Chemistry, Voi, 91, No. 19, 1987 4875
T
I 0
2
4
6
8
1
0
LAYER
-20
-15
-10
-5
0
5
1'5
1.0
2'0
Figure 2. Average number of nearest-neighbors(solid line) and hydrogen-bonded neighbors (dashed line) as a function of horizontal layer. Error bars are 1 2 estimated standard deviations.
i
I
\
-I
0
2
'
I
4
'
I
6
"
'
8
I
1
'
0
LAYER
-20
-15
-10
-5
0
5
10
15
20
(4 Figure 1. Average density profiles for hydrogen, pH(z), and oxygen, po(z), atoms relative to their bulk densities in liquid water at 1.0 g/cmS. Error bars are f 2 estimated standard deviations.
large systems is reasonable. However, in order to achieve a well-controlled description of aqueous solution surfaces, it is important to test the accuracy of simpler treatments by performing directly comparable calculations which do not truncate the Coulomb interactions. Furthermore, it is of interest to inquire how dielectric behavior is correlated with surface properties.20~2' Thus, it would be helpful to have Coulomb forces evaluated by Ewald summations in the subsequent work on these systems.
Results The atomic density profiles are shown in Figure 1. The density of the water in the middle of the lamella is roughly 0.97 g/cm3. Within the indicated statistical uncertainty, these profiles display no spatial layering of the 0 and H atoms. Furthermore, the differences between the oxygen and hydrogen profiles are zero to within the statistical uncertainty for the difference. (We discuss below the dipole density and the electric field profile which can be related to the differences in the 0 and H atom densities.) In contrast, the 0 and H atom density profiles obtained from calculations on water next to a planar wall (even a wall not wetted by the liquid) do show some significant differences and do indicate some atomic layering.* The subsequent analysis of molecular structure as a function of altitude with respect to the lamella used a spatial division of the system into layers 1.5 A in depth. The total depth of the lamella accommodated 20 such layers with appreciable average molecular population. Since the average properties of the system were accurately symmetric about the midplane (z = O.O), these 20 layers were combined into symmetry-equivalent groups of two in presenting the results below. These layers were sequentially numbered with the outermost layer designated layer 1 and the innermost layer labeled layer 10. Note that the average molecular population of layer 1 was only 0.41 molecules. Thus the average properties of molecules in this outermost layer at statistically less ~
(20) Sullivan, D.E. J . Chem. Phys. 1975,63,3684. (21) Fulton, R.L.J . Chem. Phys. 1976.64,1857.
Figure 3. Average normalized excess z component of the dipole moment of each horizontal layer normalized by the molecular dipole moment. Error bars are 1 2 estimated standard deviations.
certain than those of molecules in other layers. The average molecular population was essentially the same in layers 5 through 10, about 23 molecules. The interfacial region can be considered to be spanned by layers 1 through 5. Figure 2 shows the average number of near-neighbors and average number of hydrogen-bonded neighbors of molecules in each layer. The number of near-neighbors of a particular molecule was defined as the number of other molecules with an oxygen atom within 3.2 A of the oxygen atom of this molecule. The number of hydrogen-bonded neighbors of a given molecule was calculated as the number of near-neighbor molecules which contributed a pair energy with this molcule lower than -2.0 kcal/mol. The average number of near-neighbors is 3.9 in the middle of the lamella. This value decreases through the interfacial region to about 1.6 in layer 1. The number of hydrogen-bonded neighbors behaves similarly, falling from a value of 3.2 in the inner layers to 1.6 in the outermost layer. All near-neighbors in the outermost layer are also hydrogen-bonded neighbors. Figure 3 presents the average net z component of the dipole moment of the molecules in each layer. The net dipole moment of each layer was normalized by the magnitude of the dipole moment of the water molecule. The positive z direction was taken to be from the exterior (vapor) to the interior (liquid) of the lamella. The molecules in the outermost layers 1 and 2 do not experience a significant dipolar ordering. Layers 3 through 5 have an average net z dipole moment of between 1 and 2 multiples of the molecular dipole moment. These layers are on the inner side of the interfacial region and have average molecular populations of between 11 and 22 molecules. Thus the observed average net dipole of these layers could be regarded as a 5-10% bias of an isotropic orientational distribution of single-molecule orientations. The average dipole density in the layer with maximum polarization was 1.4 X lo-* C/m2. The excess net dipole moment points toward the interior of the lamella. This result agrees qualitatively with the expectation of Ben-Naim and Stillinger.' An independent assessment of the accuracy of our calculation is given by the fact that the x and y components of the average dipole moment in every layer is zero to within the estimated statistical uncertainty. For each interfacial layer we evaluated the probability, p ( 8 ) , of the angle, 8,formed between the molecular dipole and the
4876 The Journal of Physical Chemistry, Vol, 91, No. 19, 1987 r
'
i
i
i
i
Wilson et al.
i
0.010
om/ LAYER 4
i
0 LAYER3
t
0,005
E
LAYER 3
o1
LAYER2
0.010
0.02
t
T
-
LAYER1
4
1
0.005 I
Ib
0
0.015
0
30
60
90
120
150 180
0 (degrees) Figure 4. Probability distribution of the angle between the molecular dipole and the inward-pointingsurface normal ( z axis) as a function of horizontal layer (see eq 1). Error bars are f2 estimated standard de-
viations. inward pointing surface normal. These probabilities are shown in Figure 4. The polar angular interval from 0' to 180' was divided into 12 segments of equal length, A 8 = 15'. The average number of molecules with a dipolar angle lying within the mth angular segment was denoted by (n,,,). The probability at the midpoints, e, of the angular segments, 18 - 0,l I A8/2, was estimated from the relation [The normalization actually adopted for these figures is A8Cmp(8,) = 1.] These distributions show no striking molecular structure but do indicate that likely orientations of a water molecule in the interfacial region are with dipole lying roughly parallel to the interface. The net z component of the dipole moment in these layers is reflected in a slight positive excess of the population on the left side of these graphs (near 8 = 0') over that on the right side (near 8 = 180'). Figure 5 shows p ( 0 ) obtained for angles 8 which OH bonds make with the z axis. Since each molecule has two OH bonds which form an angle of 104.S0, it is natural to anticipate a bimodal structure for this p ( 0 ) with distinct maxima separated by an angle of 104.5' or less. Such a structure was seen by Lee, McCammon, and Rossky*for water near a hydrophobic wall. A similar bimodal distribution is observed in the present results for layer 2. The results for layer 1 are consistent with such a structure, but the large statistical uncertainties make the precise situation unclear. An interpretation of the bimodal structure of p ( 8 ) for layer 2 is that a probable orientation for molecules in this layer has one OH bond pointing outward (toward the vapor) and the second OH bond making an acute angle of roughly 70' with the z axis. Note, however, that in layer 2 the fraction of OH bonds making an angle of 120' or less is about 80%. Therefore, an appreciable fraction of the time a hydrogen bond accepting site (or "lone pair") must be oriented outward. In this circumstance, two O H bonds would form an angle of approximately 70' with the inward-pointing surface normal. The presence in layer 2 of both of the structures identified above with roughly equal probabilities is consistent with our earlier result that the average dipole moment in layer 2 is small. However, it should be borne in mind that this interpretation of p ( 8 ) is based on an idealized ordered structure of the water interface, while the real p ( B ) is quite diffuse. For layers 3 through 5 no bimodal structure was observed. A bimodal pattern which was inverted in inner layers relative to that of the outer layers (as, for example, in Figure 5 layer 2) was seen by Lee, McCam-
0
30
60
90
120 150 180
0 (degrees)
Figure 5. Probability distribution of the angle between OH bonds and the inward-pointingsurface normal ( z axis) as a function of horizontal layer. Error bars are f 2 estimated standard deviations.
mon, and Rossky.* When such an inverted pattern is observed, it is natural8 to interpret it by appealing to the variation of orientational order along the crystallographic c axis of hexagonal ice 1: if one layer has an OH bond parallel to the c axis and pointing outward, then a nearby layer can have an OH bond parallel to the c axis and pointing inward. In contrast to this interpretation, the present results for the inner interfacial layers suggest a mild preference for the OH bonds to orient parallel to the interface with a slight statistical preference for orientations with OH bonds pointed inward over orientations with the OH bonds pointed outward. The difference of the electrostatic field of a water molecule from a dipolar field contributes to the differences between the distribution of angles of molecular dipoles and molecular OH bonds. One way to focus on these differences is to calculate the mean electrostatic field at particular positions in the interfacial region. The x and y components of this mean electric field E are zero in the present lamellar geometry, and E, depends only on the altitude, z ; thus E = (O,O,E,(z)). The line integral of the electric field through the interfacial region yields the electric potential difference between the two phases +(z")
- +(z? = -
~ ~ ' k , ( zdz)
(2)
a quantity of general interest. A method for estimating the desired electric field can be established by recalling that the electric field due to a uniform charge distribution on a plane z = zo with surface density q = Q / A is
E,(z) = 2aq sgn ( z - zo)
(3)
Therefore, given a charge distribution which depends only on z , the electric field at a particular plane can be obtained by computing the difference between the total charge per unit area, respectively, below and above the plane in question. Consequently, we estimated the electric field from the formula ( E A z ) ) = 2 ~ ( Q - ( z-) Q + ( z ) ) / A
(e-(.))
(4)
where Q+(z) is the total charge above (below) the xy plane at altitude z. The x y cross-sectional area of the simulation cell is A . Since water molecules are electrically neutral units, a particular configuration makes a nonzero contribution to ( E z ( z ) only if the plane cuts the electrical charge distribution of individual
The Journal of Physical Chemistry, Vol. 91, No. 19, 1987 4877
Molecular Dynamics of the Water Liquid-Vapor Interface
L
4
-I
..
;'
-I
/
t I
1
'-/ 1
-16
I
I
I
-12
I
-8
I
I
I
-4
z(A) Figure 6. Electric fields in the water liquid-vapor interface. The results have been folded about the symmetry plane z = 0. The points are the electric fields estimated according to eq 4. The estimated standard deviation at each point is typically about half the maximum mean electric field. The solid line is the data after filtering as discussed in the text. The dashed line is the electric potential derived from the smoothed curve with the potential of the vapor phase set equal to zero. The scale gradations are gigavolt/meter for the electric field, and decivolt for the electric potential.
molecules. In the vapor phase this happens almost never with our data, and the electric field is small. Furthermore, in any isotropic bulk phase the average charge partitioned above and below the plane is equal, and the average electric field is zero. In eq 4 we used the partial charges of the TIP4P which were chosen to provide a reasonable description of intermolecular forces. Our results on the electric fields in the model aqueous film are shown in Figure 6. The data are very noisy. The estimated standard deviation of the data is typically about half of the maximum mean electric field. The magnitude of the maximum mean electric represents a charge imbalance (lQ+(z)I) of roughly one-half the charge of the electron over the area (21.71 A)2. The natural coarse-graining implicit in the usual measurement of fields, for example, the binning involved in the density profiles of Figure 1, is not present in the formula by which the mean electric field is estimated, eq 4. Therefore, we have performed a smoothing of the estimated mean electric field by a spectral filtering method.22 A smoothed result is drawn as the solid line in Figure 6 in order to guide the eye, and in an attempt to emphasize the cruder longer wavelength features which are likely to have greater statistical significance than the short wavelength features. The electric potential which may derived from the smoothed electric field results is shown as the dashed curve in Figure 6. The smoothed electric field result and the unsmoothed data indicate an average electric field in the interfacial region which points from the vapor to the liquid. This electric field is parallel to the polarization which was found above. (See Figure 3.) Despite the fact that the dipole density points inward, a test charge approaching the interface from the vapor side encounters and passes a positive charge layer first. This indicates that on the length scales established by the density inhomogeneity, the electric field sources must be considered as significantly nondipolar. For example, a molecular orientation which points the molecular dipole moment and one OH bond toward the liquid but the other OH bond toward the vapor would contribute positively to a parallel average dipole density and electric field. This issue can be considered from a macroscopic point of view also. Over large enough length scales that only dipolar charge distributions must be considered we expect that D = E 47rP.
+
(22) Matsouka, 0.;Clementi, E.; Yoshmine, M. J . Chem. Phys. 1976.64, 1351. Lie, G. C.; Clementi, E.; Yoshmine, M. J . Chem. Phys. 1976.64, 2314.
In our case the displacement field is zero, D = 0, so that E = -4aP, antiparallel to the polarization. This is not what is found, so we conclude that the charge distribution varies too rapidly in the interface region to limit consideration to just the dipole density
P. On the inner (liquid) side of the interface the smoothed result describes an electric field which changes sign and then gradually tails toward zero. The unsmoothed data are scattered about the smoothed curve; however, it should be noted that the unsmoothed data do suggest a molecular scale structural oscillation of the mean electric field in this inner region. This structuring does not survive the smoothing shown, and statistical precision is not sufficient to establish whether this molecular scale structuring is correct. The electric potential which may be derived from the smoothed electric field result is plotted as the dashed line in Figure 6. It shows only a limited region of flat behavior in the interior of the lamella. This suggests that subsequent works on water interfaces which particularly investigate these electrical properties should consider treating thicker films. The smoothed electric potential also shows a minimum just inside the interfacial region. This reflects the fact that the smoothed electric field result passes through zero just inside the interfacial region. On a molecular thermal scale these changes in electric potential are very significant, amounting to several times kBT for a charge comparable in magnitude to that of the electron. Therefore, the statistical significance of these behaviors deserves verification by more extensive calculations. The intermolecular potentials in common use for computer experimental studies of aqueous solutions, with the possible ex~ ! ~ ~been empirically adjusted ception of the MCY p ~ t e n t i a l ? have to provide reasonable accounts of a subset of properties of bulk water. Therefore, it is a nontrivial test of these potentials to compute surface properties which can be directly compared with measured quantities. One such property is the surface tension.25 The surface tension for TIP4P water at 325 K can be calculated from our data in the following way: The contribution of the intermolecular forces to the average stress, An,", is obtained from Anau=-V1(
E
m=molecules
r,(m) F A m ) )
(5)
Here Vis the volume of the container, r,(m) is the cuth component of the position of the oxygen atom of molecule m , and Fu(m)is the uth component of the total force on molecule m. The sum is over all molecules. (The kinetic energy or ideal part of the equilibrium stress is isotropic and does not contribute to the surface tension. See eq 6 below.) This formula for the excess stress applies to rigid molecular units and can be derived by imagining that the containing vessel rigidly restricts the position of the oxygen atoms only. A calculation of the Helmholtz free energy change to first order with respect to a strain of the confining vessel leads to the required expression for the stress. The surface tension, y,can then be obtained from
where again A is the area of the surface. This formula accounts for the fact that the lamella has two surfaces. The present calculation with the TIP4P potential produces a value of (132 f 46) dyn/cm at 325 K. The quoted uncertainty is f twice the estimated standard deviation. The value for experimental water at 325 K is 68 dyn/cmeZ6 (23) Method of smoothing the data was that given in Section 13.9 of Press, W. H.; Flannery, B. P.; Teukolsky, S. A,; Vetterling W. T.Numerical Recipes; Cambridge University Press: New York, 1986. (24) Impey, R. W.; Madden, P. A.; McDonald, I. R. Mol. Phys. 1982,46, 513. (25) Reference 4 reports the value of 97 dyn/cm for ST2 water at 298 K. However, the quantitative significance of this result is unclear because eq 17 and 18 of ref 4 are incorrect adaptations of the method proposed in: Miyazaki, J.; Barker, J. A.; Pound, G. M. J . Chem. Phys. 1976,64, 3364. Note further that the latter method generally requires equation of state information to correct the computed free energy difference to the desired surface free energy (Wilson, M. A,, unpublished results). We thank H. L. Scott and J. A. Barker for discussions of these points.
4878
J. Phys. Chem. 1987, 91. 4878-4881
Conclusion This paper has presented results of molecular dynamics calculations on the equilibrium interface between liquid water and its vapor at 325 K. For the TIP4P model employed the average surface dipole density points from the vapor to the liquid. This agrees with the prediction of Stillinger and Ben-Naim. However, we expect this property to be sensitive to the detailed form of the intermolecular potential adopted for the calculation. Indeed, differences among the previous theoretical treatments of this property include choice of different sign for the axial quadrupole moment of the water molecule (referenced to a center chosen particularly for the purposes of those previous works). It appears that the current nonlinear optical methods of second harmonic generation at surfaces should be able to experimentally probe the orientational ordering of water molecules near a water liquid-vapor interface.27 However, available experimental results indicate that the observed effects are likely to be Under the conditions of the present calculation the width of the interfacial region which supports significant dipole density is only slightly larger than one molecular spacing.30 The orientational ordering of water molecules at the liquidvapor interface has been analyzed in detail. Likely orientations of water molecules have the C,, molecular axis roughly parallel
to the interface. However, the distributions are quite broad, and it is possible to interpret these results in terms of broadened intermolecular structures found in bulk water phases. It might be supposed that the breadth of these distributions is due to breaking of hydrogen bonds. However, since all near-neighbor pairs in the outermost interfacial layers are also hydrogen-bonded pairs (according to the criterion adopted), we conclude instead that the flexibility of the free surface provides a major contribution to the breadth of these distributions. The mean electric field in the interfacial region is parallel to the polarization. This indicates that the charge distribution varies too rapidly in the interface region to limit consideration to the dipole density P only in macroscopic discussions of electrical properties of this interface. The computed surface tension, 132 f 46 dyn/cm, is subject to substantial statistical uncertainties. However, this precision is sufficient to indicate that the TIP4P potential which is satisfactory for simulation of bulk liquid water does not yield surface thermodynamic properties in close agreement with experimental results for surfaces. Good tests of the suitability of intermolecular potentials for the theoretical study of surfaces, and effective modifications of available water potentials to bring about better agreement with measured surface properties, would be helpful directions for further research.
(26) Jasper, J. J. J . Phys. Chem. Ref. Data 1972, I , 841. (27) Shen, Y. R. J . Vac. Sci. Technol. 1985, B3, 1464. (28) Rasing, T.; Shen, Y.R.; Kim, M. W.; Valint, Jr., P.; Bock, J. Phys. Rev. A 1985. 537. 3 1. (29) Hicks, J. M.; Kemnitz, K.; Eisenthal, K. B.; Heinz, T. F. J. Phys. Chem. 1986, 90, 560. (30) It is now widely accepted that widths of fluid interfaces should grow indefinitely as the surface area tends to infinity in the absence of gravity- See ref 16 or: Rowlinson, J. S.; Widom B. Molecular Theory of Capillarity; Clarendon: Oxford, 1982. Therefore the conditions of computer simulations such as this one must be carefully noted in interpreting the interfacial widths found.
Acknowledgment. We thank Dr. R. D. MacElroy for his support of these calculations. This work was supported by NASA-University Consortium Interchange No. NCA2- 108, and a grant from the Research Laboratories of the Eastman Kodak Co. The computer facilities used in this work were provided by the Numerical Aerodynamic Simulation (NAS) program. The analyses of electric fields was motivated by questions raised by Frank H. Stillinger regarding the electric potential difference between water liquid and vapor phases. Registry No. HzO, 7732-18-5.
Solvent Viscosity Effects on the Rate of Side-Chain Rotational Isomerization in a Protein Molecule Indira Ghosh and J. Andrew McCammon* Department of Chemistry, University of Houston, Houston, Texas 77004 (Received: November 14, 1986)
The activated rotation of the tyrosine 35 ring in bovine pancreatic trypsin inhibitor has been simulated in model solvents that have extremely different viscosities but that are otherwise identical. Both simulations are at 300 K, but one solvent corresponds to liquid water and the other to a hypothetical glassy water. Although the ring is located in the surface region of the protein, the “freezing”of the solvent reduces the rate constant for rotation by only 50%. The time required to complete individual transitions is somewhat lengthened, apparently due to slowed relaxation both of ring-solvent interactions and of the conformation of the protein matrix that surrounds most of the ring. The slowing of these relaxations also leads to the reduction in the rate constant: the persistence of ring environments that favor rotation increases the likelihood that the ring will return to the transition-state region soon after passing through it, rather than being trapped in a stable state.
Introduction Following the pioneering work of Rahman and Stillinger, molecular dynamics simulation studies have provided a wealth of information on the microscopic structure and dynamics of liquid water and aqueous solutions.’-10 One recent line of development
in this area is the study of large biological molecules in water. Molecular dynamics simulations are now providing information that ranges from dynamic characteristics of water molecules at the surfaces of proteins to the thermodynamic contributions of solvation in the binding of ligands to genetically engineered enzymes.” Such studies are deepening our understanding of the
( I ) Rahman, A.; Stillinger, F. H. J . Chem. Phys. 1971, 55, 3336. (2) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1972, 57, 1281. (3) Rahman, A.; Stillinger, F. H. J . Am. Chem. SOC.1973, 95, 7943. (4) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1974, 60, 1545. (5) Rahman, A.; Stillinger, F. H . Phys. Rev. 1974, AlO, 368. (6) Rahman, A. In Report of the Molecular Dynamics Workshop; CECAM: Orsay, 1974.
(7) Wood, D. W. In Water: A Comprehensive Treatise, Vol. 6, Franks, F., Ed.; Plenum: New York, 1979; p 279. (8) Jorgensen, W. L. J . Phys. Chem. 1983, 87, 5304. (9) Rossky, P. J. Annu. Reo.. Phys. Chem. 1985, 36, 321. (10) Pratt, L. R. Annu. Reo. Phys. Chem. 1985, 36, 433. (1 1) McCammon, J. A,; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, 1987.
0022-3654/87/2091-4878$01.50/0
0 1987 American Chemical Society
