J. Phys. Chem. 1991,95,6211-6217
6211
4. Summary
Figure 11. Forces (ab initio) in formic acid due to an external point charge of 0.1 e, 5 A from the oxygen.
negative at the expense of the third heavy atom (through its ?r bond to the carbon), and at the expense of H, (through the u interaction). Likewise, elongating the single bond between the carbon atom and the heavy atom (N or 0)withdraws charge from the carbonyl oxygen through the r system (see VII, VIII, and IX, X). Another similarity is observed for the CH bond stretch which in both cases polarizes the carbonyl group in the same way, namely, negative charge flows from the oxygen to the carbon, which is also similar to the behavior in formaldehyde. The 0-H stretch in formic acid has a local effect and does not induce a significant flux even along the OH coordinate itself (unlike the HzO case), and, similarly to the N-H in formamide the hydroxylic hydrogen is not significantly affected by stretching other bonds. As observed for HOH in water, opening the C-O-H angle in formic acid increases the polarity of the 0-H bond.
o=c-0 IX
-
o~c=o+ X
Figure 10 shows the forces on the polar atoms of formic acid as an external point charge of +0.1 e in the molecular plane approaches the carbonyl group. Both the ab initio forces and the forces calculated with the (q#,mjm} set of electrostatic parameters are shown. The ab initio forces at 5 A are further illustrated in Figure 1 1. The general pattern is similar to that in formamide. The force on the carbonylic oxygen is practically central (pairwise point charge interaction). The carbon atom, being at the junction between two bonds, is subject to an enhanced force (larger than the force on the oxygen, which is closer to the source). The alcoholic oxygen, like the nitrogen in formamide, experiences a very noncentral force, because it is further removed from the source and separated from it by a “structured” medium that modifies the force significantly. The effect is smaller for the alcoholic hydrogen since the charge flux associated with this nucleus is relatively small.
We have shown in this work that electrostatic forces in molecules cannot be represented as arising from simple central Coulombic interactions, as commonly assumed in current force fields. Molecules are not collections of isolated atoms, interacting through empty space. Rather, molecules are structured media that relay the electrostatic interaction in a way that parallels the electronic structure and geometrical arrangement of the nuclei. Consequently, electrostatic forces in general are not isotropic. The anisotropy has a static part due to the nonspherical charge distribution around the nuclei, which is represented by atomic dipoles or higher order multipoles, and a dynamic part which is e x p r d during the vibrational motion of the molecules. Both on theoretical ground and on the basis of the cases discussed, the dynamic anisotropy is in general expected to be larger than the static anisotropy. It was found that charge flux is the main mechanism by which the central force is modified and anisotropy results. Because of this charge flux the forces are not local and depend on the response of other sites to the external perturbation. Finally, we have shown that the FR atomic multipoles reproduce the ab initio forces very well for the cases discussed. The results obtained in this work suggest that unsaturated compounds in general show large charge fluxes. The conclusion is less general for saturated systems. The two cases discussed here, HF and H20, behave differently in that respect. However, both are highly polar and one may plausibly assume that, for saturated nonpolar groups, such as methylene groups, charge flux is small.’ The systems dealt with in this work are all planar. Similar nonplanar molecules are not expected to behave qualitatively differently. For example, the C 4 stretch that induces a charge flux in the OCN moiety in formamide will do so in acetamide and N-methylacetamide as well. C h a m and K~-imm’~ have found that the polar tensor is fairly transferable among these molecules. Thus the qualitative behavior of electrostatic forces found in this work is likely to be general in nature. Further work is required, however, before a quantitative description of electrostatic forces in nonplanar molecules becomes possible.
Acknowledgment. This research was supported by the Basic Research Foundation administered by the Israel Academy of Sciences and Humanities. The author is indebted to Dr. Marvin Waldman from Biosym Technologies for many useful discussions. Registry NO.CO,630-08-0;HF,7664-39-3; H20,7732-18-5; H2C0, 50-00-0; HCONHz, 75-12-7;HCOOH, 64-18-6.
A Flexlble/Polarlzable Slmple Polnt Charge Water Model S.-B. Zhu, S.Yao, J.-B. Zhu, Surjit Singh, and G. W. Robinson* SubPicosecond and Quantum Radiation Laboratory, Departments of Physics and Chemistry, P.O. Box 4260, Texas Tech University, Lubbock, Texas 79409 (Received: March 22, 1991)
The purpose of this paper is to report the properties of a flexible water model that is computationally simple enough SO that it can be applied to the molecular dynamics study of partial ensembles of liquid water in the presence of local perturbations, such as those arising from ions, surfaces, and electrical fields. This model is based on a 3-point version with both flexible geometry and polarization. The positions of the point charges are fixed, but their values are variable according to the local field. By use of this model, we have investigated a molecular dynamics simulation of the static and dynamic properties of liquid water at room temperature and normal density. Comparisons with experimental data and other model calculations are made.
1. Introduction ~~~i~~ the past
decades, a great number of computer simulation studies concerning the structure and properties of liquid Until recently, most of these have water have been 0022-3654/91/2095-6211$02.50/0
been carried out within the framework of pairwise additive intermolecular potentials. The influence of electrostatic induction ( I ) Barker, J. A.; Watts, R. 0. Chcm. Phys. Lett. 1969, 3, 144.
0 1991 American Chemical Society
6212 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 from neighboring molecules was taken into account in an average way by the introduction of fixed effective electrical charges discretely distributed at certain points in the molecular frame. This scheme has been a reasonable approximation for simulating pure water or aqueous solutions of nonpolar molecules, where the molecular environments are spatially isotropic, or nearly so. However, in other cases, such as water near surfaces or electrically charged p l a t e ~ , 3water ~ * ~ in external fields:' or water in salt solutions,'24 explicit inclusion of dynamic electric polarization
Zhu et al. TABLE I: Force Constantsa (mdyn/A)
GI
c12
cI 3
CZl
c 2 3
c3
9.29
0.73
-1.42
9.29
-1.42
2.23
uc11
=
c12,
C,I = CI3, c, = c,,.
TABLE 11: Parameters
Function A, kcal/mol-A'2
Used for the Intermolecular Potential
B, kcal/mol.A6
695 000 (2) Rahman, A.; Stillinger, F. H. J . Chem. Phys. 1971, 55, 3336. (3) Ben-Naim, A.; Stillinger, F. H.In Srruc~uresund Trunspori Prixesses in WcUcr und Aqunnw Soluffom; Home, R. A. Ed.;Wiley-Interscience: New York, 1972. (4) Popkie, H.;Kistenmacher, H.; Clementi, E. J . Chem. Phys. 1973,59, 1325. (5) Watts, R. 0. Mol. Phys. 1974, 28, 1069. (6) Stillinger, F. H.; Rahman, A. J . Chem. Phys. 1974, 60, 1545. (7) Lemberg, H. L.; Stillinger, F. H. J . Chem. Phys. 1975, 62, 1677. (8) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976,64, 1351. (9) Lie, G.C.; Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976,64,2314. (IO) Watts, R. 0. Chem. Phys. 1977, 26, 367. ( I I ) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1978, 68, 666. (12) Impey, R. W.; Klein, M. L.; McDonald, 1. R. J . Chem. Phys. 1981, 74, 647. (13) Berendsen, H. J. C.; Postma, J. P. M.; von Gunsteren, W. F.; Hermans, J. In InIermolecuIur Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981. (14) Jorgensen, W. L. J . Am. Chem. Soc. 1981,103,335; J . Chem. Phys. 1982, 77, 4156. (15) Reimm, J. R.; Watts, R. 0.; Klein, M. L. Chem. Phys. 1982,64,95. (16) Mom, M. D.; Rice, S.A. J . Chem. Phys. 1982, 76,650. (17) Jorgenscn, W. L.;Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J . Chem. Phys. 1983, 79, 926. (18) Bopp, P.; Jancso, G.; Heinzinger, K. Chem. Phys. Lerr. 1983,98, 129. (19) Reimers, J. R.; Watts, R. 0. Mol. Phys. 1984, 52. 357. (20) Reimers, J. R.; Watts, R. 0. Chem. Phys. 1984, 85, 83. (21) Reimers, J. R.; Watts, R. 0. Chem. Phys. 1984, 91, 201. (22) Toukan, K.; Rahman, A. Phys. Rev. 1985, 318, 2643. (23) Lybrand, T. P.; Kollman, P. A. J. Chem. Phys. 1985, 83, 2923. (24) Dcmontis, P.; Suffritti, G. B.; Fois, E. S.; Gamba, A. THEOCHEM 1985, 21, 201. (25) Neumann, M. J. Chem. Phys. 1986.85, 1567. (26) Lie, G.C.; Clementi, E. Phys. Rev. 1986, 33A, 2679. (27) Teleman, 0.; Jansson, B.; Engstrbm, S. Mol. Phys. 1987, 60, 193. (28) Costa, B. J.; Rivail, J. L.; Bigot, B. J . Chem. Phys. 1987, 86, 1467. (29) Ullo, J. J. Phys. Rev. 1987, 36A, 816. (30) Brodskaya, E. N.; Rusanov, A. I. Mol. Phys. 1987,62, 251. (31) Anderson, J.; Ullo, J. J.; Yip, S.J . Chem. Phys. 1987, 87, 1726. (32) Evans, M. W.; Lie, G.C.; Clementi. E. J . Chem. Phys. 1988, 88, 5157. (33) Sprik, M.; Klein, M. L. J . Chem. Phys. 1988.89, 7556. (34) Watanabe, K.; Klein, M. L. Chem. Phys. 1989, 131, 157. (35) Zhu, S.-8.; Robinson, G.W. Proc. fnr. Con$ Suprcompur. 1989, I f , 189. (36) Kuwajima, S.; Warshel, A. J . Phys. Chem. 1990,94, 460. (37) Ahlstram, P.; Wallqvist, A.; Engstrh, S.; J6nsson, B. Mol. Phys. 1989,64563. The model for comparison (PSPC) is an optimized polarizable,
rigid SPC model with enhanced repulsion. For notational consistency, we will
call this model SPC-P. (38) Cieplak, P.; Kollman, P.; Lybrand, T. J . Chem. Phys. 1990,92,6755. (39) Lee, C. Y.; McCammon, J. A.; Rossky, P. J. J . Chem. Phys. 1984, 80, 4448. (40) Zhu. S.-B.; Robinson, G. W. J. Chem. Phys. 1991, 94, 1403. (41) Evans, M. W.; Lie, G.C.; Clementi, E. J . Chem. Phys. 1987, 87, 6040. (42) Engstrh, S.; Jansson, B.; Impey, R. W. J. Chem. Phys. 1984,80, 5481. (43) Grim, 0. A.; Haymet, A. D. J.; Banet, M. J.; Simon, J. D. J . Phys. Chem. 1988, 92, 3391. (44) Maroncelli, M.;Fleming, G. R. J . Chem. Phys. 1988, 89, 5044. (45) Zhu, S.-B.; Lee, J.; Zhu, J.-B.; Robinson, G. W. J . Chem. Phys. 1990, 92, 5491. (46) Zhu, S.-B.; Robinson, G. W. Z. Nururforsch. A 1990, 16a 221. Zhu, S.-B.; Zhu, J.-B.; Lee, J.; Robinson, G. W. Ionic Dissociation Dynamics in a Polarizable Liquid, J . Mol. Llq., submitted.
600
4%. ea -0.65
A.
ea +0.325
a,A3 1.271
Oe is the electron charge.
effects seems necessary. In some recent publications,21~23~28~33~35-38~47 attempts have been made to incorporate many-body polarization effects into the simulations. However, lengthy iterations3' or other complications, which add much effort to the already time-consuming algorithms, are often unavoidable. Because of this, statistically meaningful studies of partial ensembles of polarizable water perturbed by solutes or surfaces have marginal tractability. If it is desired to include flexible intramolecular bonding in addition to the polarization, the computational requirements become even more demanding. In this paper we choose a sufficiently simple basic molecular dynamics (MD) model of water so that the combined effects of intramolecular bonding flexibility and instantaneously responsive electric polarization can be investigated. The basic model used is derived from the 3-site point charge model of Berendsen et which, though computationally simple, has proven to give a reasonably realistic representation of liquid water in various respects. Since intermolecular potentials can be sensitive to both intramolecular charge disturbances and geometry variations, this molecular flexibility gives rise to non-pair-additivecontributions to the intermolecular interactions. Our purpose of studying such models is not so much to continue refining them until the ultimate 3-site model is reached, but rather to investigate in a computationally tractable way the effect of such modifications on reasonably realistic water models. The inclusion of bond flexibility and dynamic polarization in the same model may provide clues as to the manner in which current water computations can be improved. An equally important goal is to create a simple but reasonably realistic model for a waterlike substance that may be practical in biodynamics calculations and in investigations of aqueous systems containing ionic or polar solutes. A description of the model is given in section 11, while section 111 presents the results and compares them with those obtained from other model calculations and experimental measurements. 11. Model
In this work we describe a series of MD calculations on a recently developed simple point charge (SPC)water model, a b breviated SPC-FPto indicate that both flexible bonds (F) and polarization (P)are included. This model is based on the simplest of 3-site version^,'^.'^ where one point represents the oxygen atom carrying a negative charge, and the other two represent the hydrogen atoms having positive charges. These atoms are linked to one another through purely harmonic potentials. The equilibrium 0-H bond length and H-0-H bond angle for an isolated unperturbed SPC-FPmolecule are taken to be the experimental gas-phase values:" boH = 0.9572A and $I = 1 0 4 . 5 2 O . One goal of our computations is to see how the intramolecular structure and vibrations behave under the influence of liquid-state perturbations. We follow here the harmonic intramolecular bonding prescription of Toukan and Rahman.= Let Ab, denote the 0-H bond stretches ( i = 1 or 2) or H-H distance change (i = 3). Using (47) Sauniere, J. C.; Lybrand, T. P.; McCammon, J. A,; Pyle, L. D.
Compur. Chcm. 1989, 13, 313.
(48) Benedict, W. S.;Gailar, N.;Plyler, E. K. J . Chem. Phys. 1956, 21, 1139.
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6213
Flexible/Polarizable Simple Point Charge Water Model this notation, we can write the intramolecular potential in the following form 3
uintra
=
c CiJAblAb]
72 IJ'I
(1)
system fixed in the molecule such that the x axis is the H-H line and the z axis is normal to the plane of the molecule, the polarization is simulated by additional charges on the ith molecule determined by CY@$') - E:')) --CYEp AqHf) = (3) rOHl cos (fl,ei) rH~H2
The values of the force constants C, are given in Table I. In order to track the intramolecular motions of the hydrogen atoms acCY(^!$) - E$i)) aE2) curately, extremely short (0.25 fs) integration time steps were used in the computation. A&$) = (4) rOH2 cos (!@i) rH1H2 Similar to the cases of the ordinary rigid 3-site SPCI3 and the early TIPS4 models, the intermolecular interactions between any A& = -&HIi) - AqJ) (5) pair of water molecules are comprised of nine Coulomb forces Here E!JJ are the local electric field components from the previous and a single Lennard-Jones force centered on the oxygen atom. step, CY is the (assumed) isotropic polarizability, the r's are the In order to keep the computational effort to a minimum, the * ~ ~ ~ ~atom-atom distances, and Bi is the bond angle. This method of simulation is carried out in a reaction field g e ~ m e t r y . ~This accounting for the polarizations is similar to the equivalent torque appears to be the most economical method for treating the method used in ref 57 for CSz. long-range Coulomb interactions and results in almost identical Almost every polarizable water model so far developed treats dielectric properties compared with methods based on an the oxygen atom as the center of the induced dipole Ewald-type lattice Thus, the intermolecular interaction and employs the experimentally determined gas-phase isotropic can be expressed as polarizability of 1.444 Some ambiguity is evident here, since uintcr. the experimentally determined polarizability is principally contributed by the deformation of electronic clouds, and it would be suspected that this macroscopically averaged polarizability may not remain suitable as a measure of the degree of molecular level polarization in the condensed phase. The polarizability and the where the subscripts I and k represent 0 or H; q denotes the polarization properties must in fact be affected by the local appropriate point charge; riQkis the distance between particle I molecular environment, which can give rise to significantly varying of the ith molecule and particle k of thejth molecule; and r, = electric fields on spatial scales comparable to the molecular size. 9.863 A, half the box dimension, defines the cutoff radius beyond How to choose the center of polarization and its magnitude then which all intermolecular interactions are truncated to zero. In becomes an ambiguous question, and the accurate realization of the above equation, cRF is the frequency-independent dielectric polarization effects in a model that utilizes a configurationof point constant of the background continuum. NeumannSohas shown charges becomes difficult. Nevertheless, for simplicity, and bethat the modification made in eq 2 is completely equivalent to cause errors are moderated for distances large compared with the the Onsagers5 reaction field model, provided the molecules are molecular dimensions, we adopt the following approximation. We electrically neutral and a center-of-mass cutoff is applied to the take the polarization center to coincide with the oxygen atom, whole molecule rather than to the individual sites. In the present where the deformation of electronic clouds is much more sigcalculation, we choose cRF to be infinity, corresponding to the nificant than it is on the hydrogen atoms. In other words, the conducting boundary condition. This is acceptable for moderately electric field strength in eqs 3-5 is calculated at the location of to highly polar substances such as water. Table I1 provides inthe 0 atom. We then choose the polarizability to be an adjustable formation on the potential parameters. The point charges & and parameter in order to reflect an effective quantity. In this simqg have been selected so that the dipole moment of an isolated ulation, the optimized effective polarizability was determined to SPC-FPmolecule equals 1.83 Debye, the experimentallyobserved be 1.271 A3, about 12% lower than the experimental gas-phase value in the vapor phase of water.56 Another of our goals will value.56 be to determine how this average dipole moment is modified by The SPC-FP water model is applied to an MD simulation of the liquid environment. a constant number-volume-energy (NVE) ensemble with conTo realize the polarization contributions in the SPC-FPmodel, ventional periodic boundary conditions. The ensemble contains additional charges are introduced on the oxygen and hydrogen 256 molecules in a 19.726-A cubical box, each trajectory of which atoms. These additional charges complement the fixed charges is solved by using a third-order predict-estimate-correct algoralready described. A very simple method is used, whereby the ithmOs9 In the present work, we solve the equations of motion additional charges are adjusted at each step in the MD compuwith a time step of 0.25 fs, while calculating the more timetation according to the local electric field from the previous step. consuming, but relatively slowly varying, intermolecular forces This procedure is justified since the electric field varies slowly in at every second step. Through the use of velocity scaling, the time. Indeed, Ahlstram et ala3'have investigated the electric field system is initially equilibrated at room temperature (298 K) with autmmlation functions during a simulation of polarizable water the normal liquid water density of 0.997 g/cm3. After reaching molecules and found that a decay of only about 1% occurs after thermal equilibrium, an additional 160000 time steps (40ps) are 10 fs. This corresponds to 40 time steps in our calculation. performed on the CRAY YMP/832 to accumulate data. During Because of the ultrashort time steps used, we found that iterthis accumulation phase, no further temperature adjustments were ations were not required in order to obtain stable solutions. In found necessary. The energy was constant to within 0.1%. addition to this, complicated tensorial computations were avoided 111. Results because we explicitly consider the Coulomb forces rather than In this section, we present results based on the model introduced interactions among dipole moment vectors. By use of a coordinate above. We then compare these results with those from other model calculationsand experimental measurements. These water models 149) Barker. J. A,: Watts. R. 0. Mal. -Phvs. -1973. 23. 189. include the rigid SPC model13 and its flexible-bond modificaisoj Niuma'nn, M. ~.ChCit.Phyi1985,82, s 6 K ' t i ~ n ~ ~with J ' intramolecular harmonic vibrational degrees of (51) Ewald, P.P. Ann. Phys. (Purls) 1921, 21, 1087. freedom (abbreviated SPC-Fl). We also take as other points of (52) Brush, S.0.;Sahlin, H. Lo;Teller, E.J. Chem. Phys. 1966,45,2102. (53) Janmne, V. M. Chem. Phys. 1974,3, 78. comparison a flexible version26of the nonempirical MCY po(54) Ladd, A. J. C. Mol. Phys. 1978, 36,463. (55) Fmhlich, H. Theory of Dfelecrrlcs;Oxford University Pres: Oxford, U.K.,1958. (57) Zhu, S.-B.;Zhu, J.4.; Robinson, G. w. Mol. Phys. 1989,68, 1321.
+-
I -
(56) Eiwnberg,, D.; Kauzmann, W. The Structure a d Properties of Water, Oxford University Press: Oxford, U.K., 1969.
(58) Stillingcr, F. H.; David, C. W. J . Chem. Phys. 1978, 69, 1473. (59) Bceman, D. J . Compur. Phys. 1976, 20, 130.
Zhu et al.
6214 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 3.5
TABLE 1Ik Geometry of the Liquid-State Water Molecule experiment (gas phase)' experiment (liquid)b SPC-FP (this work) SPC-FI' ZR-FPd MCYL-F'
( m u ) ,A 0.957
0.966 0.966 1.016 0.979 0.975
(HQH), deg
3.04
I
(1'1.
104.5 102.8 101.0 104.9 102.0 103.5
'Reference 48. bReference 63. 'Reference 27. Simulation F. Reference 35, 'Reference 26. For notational consistency throughout, a suffix F has been attached to emphasize that this model uses flexible bonds.
tential: denoted for consistency in this paper as MCYL-F, and a rather more elaborate flexible-polarizable 3-site model of our labeled ZR-FP. Some TIP4P results are also compared here. The inclusion of TIP4P in the comparisons will provide a benchmark for early rigid models, since Jorgensen et al.I7 have already carried out an extensive comparison between rigid nonpolarizable 3-, 4- and 5-site water models, including SPC,TIP4P and ST2. Since the preparation of this paper began in mid-1989, some other papers on polarizable ~ater~'*~*@ or flexibly bonded watefll*62have appeared in the literature. Where possible, some of these models will also be compared with the SPC-FP model. A recent paper on an optimized polarizable, but rigidly bonded, SPC model of water,37which we will call SPC-P, is particularly relevant to our work here. To our knowledge, the SPC-FP model described here and the earlier ZR-FP modePs from our laboratory have been the only ones to include simultaneously both polarizability and flexible bonding in a water MD simulation. In Table 111, the average geometry of water molecules obtained from the SPC-FPmodel is reported and compared with data from neutron-scattering experiment^.^^ Comparison with the results from other MD simulations, which employ flexible intramolecular bonding, is also given there. It is seen that the bond length from the present simulation is in excellent agreement with the experimental value, while the calculated bond angle is somewhat too small. We have purposely refrained from comparing the thermodynamic averages from various computer simulations for two reasons. First of all, the usually quoted experimental resultsM$s are based on assumptions that are difficult to establish indepe11dent1y.I~ Furthermore, there seems to be no commonly accepted way to compare polarizable vs nonpolarizable, nor rigid vs flexible models. For the SCP-FP model at T = 298 K, the intermolecular energy U is -9.81 kcal/mol, while the pressure P i s about 1 atm and the density p is 0.997 g/cm3. Illustrated in Figure l a - c are site-site radial distribution functions (RDFs), goo(r),goH(f)r and gHH(f). In these figures, the curves obtained from the present simulation are compared with those calculated from the SPC-FI model and the results from the X-ray data of Narten and LevyMand the more recent neutronscattering experiments of Soper and Philli~s.6~The first peak of the 0 partial structure function in the SPC-FP model is not so high as in the SPC-F1 model and is much lower than that in the original SPC model (not shown in the figure), evincing the effects of the intramolecular flexibility. Also because of this, the shallow first minimum around 3.3 A is missing, and in fact the distribution function is higher there than in other models. Both the SPC-F1 and SPC-FP models are seen to give too little 0-0 structure beyond this region. Clearly, some of the excess density near 3.3 A needs to be moved out to the neighborhood of 4.5 A, which marks the location of the closest non-H-bonded oxygen (60)Wall vist, A. Chem. Phys. Lerr. 1990,lbS, 437. (61)G u d i a , E.;Padr6, J. A. J. Phys. Chem. 1990,946049. (62)Dinur, U. J . Phys. Chem. 1990,91,5669. (63)T h i w n , W. E.;Narten, A. M. J . Chem. Phys. 1982, 77. 2656. (64)Dorscy, N. E.Proprtles of Ordlnary Water Substances; Rheinhold: New York, 1940. (65)Kell, 0. S.J . Chem. Eng. Data lWS, 20,97. (66)Narten, A. H.;Levy, H.A. J . Chem. Phys. 1971,S5, 2263. (67)Soper, A. K.; Phillips, M.G . Chem. Phys. 1986, 107,47.
1 .a
1.2 LL
n E
0.6
0.0
1.0
1.5
2.0
2.5
3.0
3.5
4.0
4.5
Distance
1.0
-
L L -
n E
-
0.5 -
I I , I I I I I I , I I I I I ( I ( I I ( I I I ( I I I I (
1.5
2.0
2.5
3.0
3.5
4.0
4.5
5.0
Distance Figure 1. Radial distribution functions of liquid water at 298 K and 0.997 g/cm3: (0)SPC-Fp; (- -) SPC-FI;(-) experimental results from ref 66;(---) experimental results from ref 67. (a, Top) goo. (b, Center) gOH. (c, Bottom), g H H -
neighbor in tetrahedrally packed water molecules.s6 The flattened second peak has been thought" to be a general problem for 3-point models, and even for some14-point models. See the recent discussion by Cieplak, Kollman, and L ~ b r a n d . ) This ~ problem, even for 3-point models, by using however, can be corre~ted?~ sharper angular dependent intermolecular interactions (see below). These results seem to suggest that to correct for the missing gm(r) structure beyond the first peak, it might be worthwhile to consider a compromise between the present model and a polarizable water model with out-of-plane charge points, as in ST2.2 It is known that conventional rigid models of this type overestimate the liquid structure, giving a first 0-0 peak that is too high and sharp.6." However, it is also interesting to note that the recently reported polarizable water model (SK-P) of Sprik and Klein)) gives a 0-0 RDFs as strongly structured as those obtained from these rigid models. The radial distribution functions go&) and gHH(f)are Similar for the SPC-FP and SPC-FI calculations, but again the presence of both polarizability and flexible bonding is seen to weaken the structure compared with that from flexible bonding alone. Compare the curves using circles (SPC-FP) with those using long dashes (SPC-F1) in Figure lb,c. In all cases, the introduction of polarization reduces the structure and gives rise to lower first
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6215
Flexible/Polarizable Simple Point Charge Water Model
3.0
2.5
I
TABLE I V T n n ~ h t i oDiffusion ~l Constant (10-5 cd/sec) experiment' 2.3 SPC-FZd 2.5 SPC-FP (this work) 7.1 SPC-P 2.0 SPCb 4.4 TIP4Pf 2.8 SPC-Fl' 6.1 MCYL-P 1.9
'Reference 69 at 298.2 K. *Reference 27. 'Reference 27, harmonic potential. dReference 3 1 , Morse potential. 'Reference 37. /Reference 25, at 293 K. #Reference 26. TABLE V WelecMc Properties
1
0.0
2.5
I
! , ? , I
L
I
3.b
I
I
I
I
I
I
I
3.5
I
I
I
I
4.0
I
I
,
I
4.5
I
,
I
I
,
5.0
,5.5I
Distance
SPC'
1.8
1
L
L
T
n
SPC-F2' TIP4Pg MCYL-Fh
-
-
1
1.0
GK 4.3' 2.65 2.61 4.2 3.14 1.45
gK 2.90 1.78 1.75 2.8 2.11 0.98
8d 54.6 67.8 82.5 53 26
ps 7.3e 6.0 11 6.32 6.1
ps 2.SCd 1.3
ra,
5.0 3.25 4.1
1.0 I
0.6 -
0.0
D
2.4 2.44 2.28 2.42 2.18 2.26
'As cited in ref 25 for T = 293 K. bMeasured70for 1.00 g/cm3 isochore at 293 K. 'See text. dReference 73. Also, see text. eReferencc 74, run b at 300 K. /Reference 31, at 300 K, Morse potential. #Reference 25, at 293 K. *Reference 26, run b at 300 K.
-
1.2 L -
/&
experiment SPC-FP (this work)
1.5
2.0
2.5
3.5
3.0
1
4.0
4.5
Distance
0.8
LL
0
1 6 .
0.6
0
200
400
600
800
1( 30
Time (fs) Figure 3. Time correlation functions for the molecular vector along the bisector of HQH bond angle: (- -) SPC-FP(-) SPC-Fl from ref 27 with harmonic potential. 0.0
L/t r I
,
,
,
I , , ,
I
,
,
,
I
,
,
,
I
,
,
,
I
,
(68) Chen, S.H.;Teixeira, J. Adu. Chrm. Phys. 1986. 61,1.
,
,
Teleman et a].*' have shown that the most important effect of the inclusion of intramolecular structural flexibility is on the transport properties, which can in general be enhanced by a factor of 2-3. This conclusion has been recently confirmed by Gugrdia and Padrb.6' However, these results are at variance with the simulations of Anderson, Ullo, and Yip.)' The latter authors found that the flexibility of the water molecule reduces the self-diffusion coefficient appreciably. This seems like an unusual result, since molecular flexibility, all other things being equal,cannot help but lead to a "softening" of the intermolecular potentials and thus to faster diffusion. Though the effect seems to lie in the wrong direction, the key to this disagreement could be the choice of intramolecular potentials. Teleman et alaz7treated the SPC water molecule as a harmonic oscillator (SPC-Fl), just as we have done in the present calculations. On the other hand, Anderson and his co-workers" considered the 0-H bond stretches to be governed by a Morse potential (SPC-F2). This could affect the coupling between the intra- and intermolecular degrees of freedom. The importance of this coupling has been discussed recently by Dinuf12 with respect to the properties of water under electrostatic forces. A hint given by the various results may also be that the large self-diffusion coefficient resulting from the SPC-FPmodel could be made more realistic by modifying the harmonic potential. Among all the models listed in Table IV,the SPC-FP molecule has the smallest geometrical size and, as seen from Figure 1, the weakest hydrogen bond structure. These properties give rise to transport processes that are too fast when compared with ex~eriment.6~The effect, caused by the %oftcning" of the in(69) Krynicki, K.; Green, C. D.; Sawyer, D.W.Faraday Discuss. Chrm. Soc. 1978,66, 199.
Zhu et al.
6216 The Journal of Physical Chemistry, Vol. 95, No. 16, 1991
stantaneous potentials when the molecules are geometrically flexible, is enhanced further when the values of the point charges are also variable. Figure 3 depicts for the SPC-FP model the time correlation function of the molecular vector along the bisector of the H-O-H angle. Thii curve decays faster than that for SPC-Fl water. Table V summarizes some d i e l e c t r i ~ and ~ ~ , dipole ~ ~ relaxat i ~ n ~ l properties -~’ of water in its liquid state at room temperature. Unlike pairwise additive models, which take the polarization into account in an average way by introducing fixed effective charges to match the large liquid-state dipole moment, the quantity p calculated from the non-pair-additive SPC-FP model sensitively depends upon the particular state of the water. Under the circumstances studied here, the average dipole moment is 2.44 D, which, according to Lie and Clementi,26matches the generally accepted experimental liquid-state value of 2.4 Debye. See also ref 56. This dipole moment also agrees well with the semiempirical polarizable model (2.5 D) of Finney and G o o d f e l l o ~as , ~ well ~ as that (2.48 D) from MCY with a three-body correction.” These results show that effects from collective polarization are of importance not only in describing aqueous ionic solutions46but also for pure water in condensed phases. Fixed charges that are assigned larger than normal values in order to emulate the solution-phase dipole moment may under certain types of perturbations give a distorted representation of liquid water.77 Also, the flexibility of intramolecular bonds affects the dipole moment in the liquid and must be properly included along with the dynamic polarization to provide an all-around realistic picture. A good discussion of the various dielectric properties of water to be discussed in the next three paragraphs can be found in Neumann’s paper2Son rigid, nonpolarizable TIP4P water. The finite-system Kirkwood GK factor N
GK =
(i= xI r / / b / l ) 2 / N
(6)
measures the equilibrium fluctuations of the collective dipole moment N
M = Cr, 111
(7)
The static dielectric constant can be calculated by using78
= 4rnGKr2/3kBT+ 1
(8) providing the boundary condition is assumed to be conducting. In q 8, n denotes the number density, and p2 is the mean-square dipole moment per water molecule in the liquid state. Using eq 6, we find that GK for the SPC-FP model equals 2.65. This gives 54.6 for the static dielectric constant to. Note that the type of boundary conditions employed has a strong and characteristic From the relationship between GK and influence on GK.25~79*80 the Kirkwood gk factor with conducting boundary condition^^^ 3 €0 GK (9) 2co + l g K €0
-
one obtains gK = 1.78. See Table V for the compilation of these parameters for some water models. (70) V e m a t ~M.; , Frank, E. U. J . fhys. Chem. Re/. 0410 1980,9,1291. (71) Hertz. H.G. In Wurer, A Compnhensiue Treatise;Franks, F., Ed.; Plenum: New York, 1973; Vol. 3, Chapter 7. (72) Pottel. R. In Wurer, A Comprehensive Treurise; Franks, F., Ed.; Plenum: New York, 1973; Vol. 3, Chapter 8. (73) Jonas. J.: DeFriea, T.;Wilber, D. J. J . Chem. fhys. 1976,65, 582. (74) Alper, H.E.; Levy, R. M.J . Chem. fhys. 1989,9/, 1242. (75) Finney, J. L.; Goodfellow, M.In Srrucrure und Dynamics: Nucleic Acids and Proteins, Clementi, E., Sarma, R. H.,Ed.;Adenine: New York, 1983. (76) Wojcik, M.; Clementi, E. J . Chem. fhys. 1986,84, 5970. (77) Beme, B. J.; Wallqvist, A. Chem. Scr. 1989, 29A, 85. (78) de Leeuw, S.W.: Perram, J. W.; Smith, E. R. froc. R. Soc. London 1980. A 373, 27. (79) Neumann, M.; Steinhauser, 0.Chem. Phys. krr. 1983, 102, 508. Pawley, G. S. Mol. f h y s . 1W. 52, (80) Neumann, M.; Steinhauser, 0.; 97.
LL 0.91
“ I
-‘.‘,.. .
\
I .
.-I \
1
0.7 0
..- -
50
1
100
150
I I
I
200
Time (fs) Figure 4. Dipole moment time correlation functions for conducting boundary conditions: (- -) collective correlation function for SPC-FP from the present work; (- -) single-particle correlation function from the present work; (-) collective correlation function for TIP4P water at 293 K.
-
Figure 4 compares the normalized ‘collective” time correlation functions of the total dipole moment of the sample Wt) =
(M(O)-Wt))/(@)
(10)
computed from the present work and from the TIP4P potential.l73 Of course, the correlation function @ ( t ) , which can be used to calculate the frequency-dependent dielectric constant of a model system, also depends upon the boundary c o n d i t i o n ~ . ~For ~~~~J~ the purpose of comparison, the analogous function describing reorientation of a single-particle dipole moment $ ( t ) is also included in the figure. Again, the great enhancement of intermolecular dynamical processes caused by the inclusion of bond and charge flexibility is evident. A comparison of the calculated relaxation times for the single-molecule dipole moment (7,) and for the collective dipole moment ( T ~ ) associated , with @ ( t ) in eq 10, with experimental measurements7’-” and with some other model simulations is also given in Table V. Some clarification of these relaxation times seems necessary. In the paper of Jonas, et al.?’ a ‘reorientation time” T # is discussed. A rough interpolation of their data at 1 bar pressure and at temperatures of 10 OC and 30 O C indicates that T # in H 2 0 should be about 2.5 ps at 298 K. Following Hertz?’ where is denoted T~ (not to be confused with the collective correlation time 7, discussed here and in the work of Neumann’), T~ is the correlation time of the mnd-order spherical harmonics, and is related to the rotational diffusion constant, 7S-I = 60,. Hertz then introduces what he also calls a “reorientation time” T , = 3rb which is the correlation time of the first-order spherical harmonics. He then discusses its relationship with the ‘macroscopic” dielectric relaxation time, usually called 7d, for which P ~ t t egives l ~ ~a value of 8.25 ps at 298 K. Some confusion arises in Potters discussion of the correlation function d(t) of the , decay molecular dipole moment vector ( r ( 0 ) - p ( t ) ) / ( p ( O ) 2 )whose is the singlaparticle ‘reorientation time, according to time” 7,. As Pottel states, and this is also discussed in refs 25 and 71, there is no generally valid theory relating d ( t ) and @ ( t ) or $ ( t ) , the time correlation function of the polarization noise. Thus there is no valid way of relating r1 and T , or T+ From are nearly the Onsager theory, Pottel finds that T, and 7d (wC) same. On the other hand, Hertz claims that the ratio T,/T, (in his notation ~ J T , ) lies between 1 and 3, very likely closer to 3 than to 1, and Neumann thinks that T D / T ~= T,/T, = gK9 which in real liquid water is indeed close to 3. Also, according to N e ~ m a n n ?in ~ the case where cRF = a, the collective dipole correlation time 7, is equal to the Debye relaxation time TD. Needless to say at this point, the plethora of notations and meanings of these rotational relaxation times in liquid water can create a certain amount of confusion. However, the T # 5 2.5 ps value of Jonas et al. seems a good choice for the experimental T , at 298 K, and gK7, = 7.3 ps seems a good choice for the experimental 7, 7D. This means that the usual Debye correlation time T D is close to but slightly smaller than the “macroscopic”
-
Flexible/Polarizable Simple Point Charge Water Model
TABLE VI: Intnmdceulrr Vibntioarl Frequencies (cm-I) experiment' SPC-FP (this work)
SPC-FIC
bend
sym str
asym str
1640 1802 1793
32006 3875 3841
34006 4080 3966
Reference 56. *These frequencies are very inaccurate because of complex liquid-state coupling, which gives a broadening of about 400 cm-I. cReference 27, harmonic potential.
dielectric relaxation time 7d. as already discussed in one of our own papers.81 In order to investigate the effect of the liquid milieu on the high-frequency internal modes, we analyzed the velocity autocorrelation functions for the oxygen and hydrogen atoms. With the help of the standard Fourier transform technique,=S6 we find that the intramolecular bending motion is centered at 1802 cm-', slightly higher than the frequency evaluated by Toukan and RahmanZ2using a flexible harmonically bound water (SPC-Fl) with parameters for internal degrees of freedom almost identical with those listed in Table I. (In ref 22, the authors also considered Morse potential coupling). As seen in Table 111, the Toukan and Rahman SPC-FI model gives a larger water molecule, both with respect to bond length and bond angle, than our SPC-FP model. In addition to this factor, the local perturbation of the polarization on the intermolecular interactions is a contributing reason why the frequencies of the stretching modes are higher than those obtained from the SPC-Fl model. As shown in Table VI, both SPC-F1 and SPC-FP frequencies are substantially higher than the experimental values.% Replacing the harmonic intramolecular or other type potential with a more realistic Morse p0tential,2~*~' of anharmonic potential,82would improve these vibrational frequencies ~onsiderably.~~ IV. Conclusions In this paper we have described a polarizable three-site flexibly bonded water model and applied this model to study the static ~
(81) Robinson, G. W.; Thistlethwaite, P. J.; Lee, J. J . Phys. Chem. 1986, 90,4224. (82) Dang, L. X.; Pettitt, B. M. J . Phys. Chem. 1987, 91, 3349. (83) Zhu, S.-B.;Singh, S.; Robinson, G.W. A N e w Flexible/Polarizable Water Model, J . Chem. Phys., in press. This is a 5-site FP model with a
specific H-bonding interaction and anharmonic intramolecular potentials prescribed by Dang and Pettitt.'2 The RDFs arc similar to those obtained for ZR-FP (Figures 2,a-c, this paper). While not affecting those properties that agreed with experimental data in the SPC-FP model, a number of significant improvements over SPC-FP have been made in other respects: vibrational frequencies, the HOH bond angle, Tc and 7,. However, the static dielectric constant (58.6) and the translational diffusion constant (1.6 X IV5 cmz/s) are both about 30% too low. The computational times for this new model are roughly 4X longer than for the SPC-FP model.
The Journal of Physical Chemistry, Vol. 95, No. 16, 1991 6217 and dynamic properties of liquid water, using the molecular dynamics method. Included in this study are the internal energy, the pressure, the sittsite radial distribution functions, the dielectric constant, the self-diffusion coefficient, dipole autocorrelation functions, and intramolecular vibrational frequencies. Comparisons with experimental results and other model calculations have been made wherever possible. It is found that this SPC-FTmodel reasonably reproduces the perturbed molecular geometry in the liquid state. However, the liquid structure, as evidenced by the 0-0 RDF, is improved only slightly, if at all, over that from the rigid SPC model. The O-H RDF from SPC-lT lies between the two experimental c u r ~ e s , ~ while 9 ~ ' the H-H RDF has too shallow a first minimum. A number of suggestions are made for improving this structure. The SPC-FP model provides a good average liquid-state dipole moment per molecule. It is also found that inclusion of polarization may strengthen effects from couplings between intramolecular and intermolecular motions, which, in general, can influence transport properties and certain electrostatic effects62by a considerable factor. The present simulation gives dielectric properties roughly comparable to experimental observations and other computational results. Overall, this simple polarizable model seems to describe many of the properties of liquid water within acceptable accuracy at one particular temperature. For this reason it could be a useful tool for studying a waterlike substance in strong electrical fields, in salt solutions or in biological systems, where the influence of external perturbations overshadows some of the more subtle details of liquid-state water-water interactions. The other goal of these studies was to learn how to improve existing water models with the addition of flexible bonds and polarizability. From the lessons learned here, it should be possible to create a still more realistic water model that is manageable under various perturbations without excessive computational effort.83 Such a model must hold its validity over a broad range of temperatures and thermodynamic conditions. The representation of polarization may have to be improved, and the structural flexibility has to be less stiff than that used here. These changes may help bring the transport properties into better agreement with experiment. In order to retrieve realistic radial distribution functions, a sharper angular dependent intermolecular interaction potential, as was used in ref 35, would seem to be essential.
Acknowledgment. Financial support at the SPQR Laboratory has been shared by the Robert A. Welch Foundation (D-0005, 33% and D-1094, 5%), the National Science Foundation (CHE8611381,44%), and the State of Texas Advanced Research Program (1306, 18%). Computer time was furnished by the Pittsburgh Supercomputing Center. Registry No. Water, 7732-18-5.
