Water Simulation Model with Explicit Three-Molecule Interactions

J. Phys. Chem. B 2008, 112, 8311–8318

8311

Water Simulation Model with Explicit Three-Molecule Interactions R. Kumar and J. L. Skinner* Theoretical Chemistry Institute and Department of Chemistry, UniVersity of Wisconsin, Madison, Wisconsin 53706 ReceiVed: January 31, 2008; ReVised Manuscript ReceiVed: April 8, 2008

Much effort has been directed at developing models for the computer simulation of liquid water. The simplest models involve effective two-molecule interactions, parametrized from experiment, for use in classical molecular dynamics simulations. These models have been very successful in describing the structure and dynamics of liquid water at room temperature and one atmosphere pressure. A completely successful model, however, should be robust enough to describe the properties of liquid water at other thermodynamic points, water’s complicated phase diagram, heterogeneous situations like the liquid/vapor interface, ionic, and other aqueous solutions, and confined and biological water. In this paper, we develop a new classical simulation model with explicit three-molecule interactions. These interactions presumably make the model more robust in the senses described above, and since they are short-ranged, the model is efficient to simulate. The model is formulated as a perturbation from a classical two-molecule interaction model, where the forms of the correction to the two-molecule term and the three-molecule terms result from electronic structure calculations on dimers and trimers. The magnitudes of these perturbations, however, are determined empirically. The resulting model improves upon the well-known two-molecule interaction models for both static and dynamic properties. I. Introduction An understanding of the properties of water is important for many fields of science and engineering. As such much effort has been directed at developing models for the computer simulation of water.1 The simplest models involve effective twomolecule interactions,2–13 parametrized from experiment, for use in classical molecular dynamics simulations. These models have found widespread use and, indeed, have been very successful in describing the structure and dynamics of liquid water at room temperature and 1 atm pressure. Even at this state point, however, these models underestimate the time scales for molecular rotation and hydrogen bond rearrangement.14 Moreover, a successful model should ideally be robust enough to describe the properties of liquid water at other thermodynamic points, not to mention water’s rich and complicated phase diagram (critical point, several forms of ice, etc.),15–18 water in heterogeneous situations like the liquid/vapor interface,19–30 the properties of ionic and other aqueous solutions,31,32 and those of confined and biological water.33–37 In this regard, the classical effective two-molecule interaction models are generally inadequate. These models are insufficient for two reasons. First, they are classical, while at this point it is clear that nuclear quantum effects associated with the motion of the light hydrogen atoms are quite important.38–44 One can capture some of the effects of these quantum fluctuations with an effective classical potential, as in the models described above, but then it is not necessarily reasonable to suppose that the model is robust (in the senses described above). Second, the Born-Oppenheimer potential surface for interacting water molecules is simply not pairwise additive. In other words, there are complicated many-molecule electronic effects.45–49 One way to deal with these many-molecule electronic effects is by incorporating electronic polarizability, which can be * To whom corresponence should be addressed. E-mail: skinner@ chem.wisc.edu.

implemented in many ways (point molecular or atomic polarizabilities, fluctuating charges and/or dipoles, Drude oscillators, empirical valence bond states, etc.).50–66 In these methods, then, the charge distribution of a given molecule can respond to its local environment, thereby modifying its interactions with other molecules. These polarizable models have been successful in improving the robustness of classical pairwise interaction models. Drawbacks of these approaches, however, are 2-fold. First, they are sometimes time-consuming to implement (since they usually either are iterative or involve auxiliary adiabatic variables). Second, these polarizable models do not often fully capture the complicated electronic structure effects (such as exchange and charge transfer) that occur with interacting water molecules.67 Efforts to improve this situation have been in two general directions. First, there are the ground-breaking and powerful ab initio molecular dynamics (MD) approaches (Car-Parrinello MD and Born-Oppenheimer MD), where the Born-Oppenheimer potential for a given nuclear configuration is computed, using electronic structure methods, on the fly.68–82 In principle, this approach completely solves the aforementioned electronic structure problem. In practice, however, there remain problems associated with the exchange functional, basis sets, pseudopotentials, convergence issues, and sample size (since the method is very time consuming to implement). The second general direction involves developing and parametrizing potentials with two- and three-molecule interactions from quantum chemistry calculations and from experiment.83–113 Such an approach implicitly assumes that four- and higher-molecule interactions are negligible, which appears to be the case for water.46–48 In principle, this kind of approach again completely solves the electronic structure problem, although there are still the usual issues associated with the level of electronic structure theory, basis set superposition error, and the complete basis set limit. In terms of nuclear quantum effects, it does not really make sense to add quantum fluctuations to an effective classical model

10.1021/jp8009468 CCC: $40.75  2008 American Chemical Society Published on Web 06/21/2008

8312 J. Phys. Chem. B, Vol. 112, No. 28, 2008

Kumar and Skinner

that has been parametrized from experiment (and so in that sense already has these effects incorporated into the model). On the other hand, when pursuing one of the ab initio methods that in principle produces the true Born-Oppenheimer potential, in order to do the problem correctly, it is essential to include quantum effects, and indeed path integral and centroid MD calculations have been performed on these models, leading to quite impressive results.38–44 These quantum calculations have their own limitations in terms of accuracy, computer time, and system size. One last approach to mention is the coarse-grained method.114,115 In this approach, water molecules are replaced by one or two sites, and the appropriate effective potential is obtained by force matching from an atomistic simulation (possibly using ab initio methods). Such an approach can be essential for very largescale simulations. From a theoretical perspective, the fully ab initio approaches together with quantum dynamics are clearly the most attractive and, in many ways, are the approaches of the future, but at this point, they are not feasible for large-scale simulations for complicated problems of current interest. On the other hand, while for very large problems the coarse-grained approach may indeed represent the best alternative, for other smaller problems where more detail is required, they may not be appropriate. We believe, then, there is still the need for improved atomistic classical simulation models that are robust enough to describe water in its many phases and in heterogeneous environments; such models presumably would find widespread use in biological, energy, and materials MD simulations. For such an effective classical model to be robust, it should capture, as faithfully as possible, the true quantum chemical effects of interacting water molecules. Our approach, described below, is an ab initio based but empirically parametrized effective classical model with explicit three-molecule interactions. Since these three-molecule interactions are short-ranged, the model is relatively easy to implement and quick to simulate. The resulting model improves upon well-known two-molecule interaction models for both static and dynamic properties. II. Development of New Potential Model In what follows we develop a classical simulation model with explicit two- and three-molecule interactions, whose form is determined by electronic structure calculations. As discussed above, the true Born-Oppenheimer potential is not the appropriate one for classical dynamics; rather for this, one needs an effective potential. Thus the magnitudes of the two- and threemolecule terms are determined empirically by comparing simulation results to experimental results for liquid water. To begin, we write the total interaction energy of a collection of rigid molecules as a sum of two-molecule (Eij) and threemolecule (∆Eijk) terms:49

U)

∑ ∑ Eij + ∑ ∑ ∑ ∆Eijk i

j>i

i

(1)

j>i k>j

where

∆Eijk ) Eijk - (Eij + Ejk + Eik)

(2)

Thus Eij and Eijk are the dimer and trimer interaction energies, respectively. We have assumed, as discussed above, that fourand higher-molecule terms are not important.46–48 A.. Two-Molecule Interactions. Rather than starting from scratch to develop a form for the dimer interaction energy, Eij, it is more convenient to write this energy as a reference term, Erij plus a correction ∆Eij:

Eij ) Eijr + ∆Eij

(3)

In fact, we choose the reference dimer interaction energy to be that for the TIP4P model,3 which incorporates the experimental gas-phase geometry, rendering it more appropriate for describing low-density systems, and for dimer and trimer ab initio calculations. With this choice of reference potential it is clear that there must be a correction, since the reference potential is really an effective two-molecule interaction and already implicitly includes the effects of three-molecule interactions. As we will see below, this choice of reference will allow the correction term to have a simple form. Choosing such a reference has two other advantages. First, many classical pairwise models provide a realistic description of liquid water at room temperature and therefore provide a good starting point for the development of this new model, and in that sense the new model can be considered a perturbation from this reference. Second, we use the reference model to generate representative dimer and trimer configurations for the ab initio calculations. Thus, we determine the form and parameters of the correction ∆Eij by fitting to ab initio dimer interaction energies. We first perform a standard MD simulation of the TIP4P water at room temperature and the experimental density (see below). Dimer configurations are selected randomly from all pairs of molecules with O-O distances of less than 5 Å. For each of five hundred pairs, MP2 calculations with the aug-ccpVTZ basis set116 were performed to obtain the ab initio interaction energies Ea12 (arbitrarily labeling the two molecules as 1 and 2), which ranged from -20 to +12 kJ/mol, and many of these energies are near zero. Considering only those dimers whose oxygen distances are less than 3.3 Å (the first minimum in the O-O radial distribution function for TIP4P water), that is, the more strongly interacting ones, the average of these interaction energies is -9.51 kJ/mol. On the other hand, the average TIP4P energy for the same set of dimers is -12.50 kJ/mol. Thus the average TIP4P energy is substantially lower than that from the ab initio calculations, since, as mentioned earlier, the TIP4P energies incorporate (mostly favorable) threemolecule interactions. We find that a positive correction term of the simple form

∆E12 ) ae-br12

(4)

provided a good fit (through eq 3) to the ab initio results, with a ) 1.925 × 106 kJ/mol and b ) 4.67 Å -1. r12 is the O-O distance between the two molecules. Thus the correction is an exponential repulsion. A scatter plot, for all five hundred molecules with O-O distances less than 5 Å, of the fit energies E12 to the ab initio energies Ea12 is shown in Figure 1. The rms deviation is 1.6 kJ/mol. The ab initio calculations described above do not correct for basis set superposition error and are not extrapolated to the complete basis set limit.99,117–120 Interestingly, for water dimers and trimers, with the MP2 level of theory and the aug-cc-pVTZ basis set, the uncorrected interaction energies are closer to those in the complete basis set limit than are to those corrected for basis set superposition error by using the counterpoise method.120,121 Thus for our particular problem, with our method and basis set, correcting for basis set superposition error is not warranted (and in addition is quite time consuming). One can still assess how far off our results are from the complete basis set limit. For half of the five hundred dimers selected above, we extrapolated to the complete basis set limit,120 and for those that have O-O distances less than 3.3 Å, the average interaction energy is -9.03 kJ/mol, which is about 0.5 kJ/mol higher than

Water Simulation Model

J. Phys. Chem. B, Vol. 112, No. 28, 2008 8313

Figure 1. Scatter plot of the fit energies E12 versus the ab initio energies a E12 . The diagonal line is a guide to the eye. The root-mean-square deviation is 1.6 kJ/mol.

the uncorrected results reported above. This difference of 0.5 kJ/mol can be compared to the average two-molecule correction of 3.24 kJ/mol from above. We are not particularly concerned if our ab initio results have a small error, since the magnitude of the two-molecule correction, a in eq 4, will be multiplied by a parameter λ1 to be determined empirically (see below). But we do want the exponential form and the value of b in eq 4 to be correct. To that end, we fit the complete-basis-set limit energies to the same functional form, and obtained a value for b of 4.44 Å -1, with about the same rms deviation as before. This value of b differs from the value for the uncorrected energies by about 5%. This small difference is probably not something to worry about, and in any case for the trimers, it would be very computationally demanding to extrapolate to the complete basis set limit, so we are content with proceding with the uncorrected energies for both dimers and trimers. B. Three-Molecule Interactions. Determination of the threemolecule terms proceeds in a similar fashion. We extract random trimers from the TIP4P simulation, with the requirement that at least two of the three O-O distances are less than 3.3 Å (the first minimum in the TIP4P O-O radial distribution function). Note that this is shorter than the cutoff used for dimer selection, consistent with the expectation that the three-molecule interactions are shorter ranged. For each trimer, we perform MP2 calculations of Eijk and Eij and Eik and Ejk so as to obtain the three-molecule energies ∆Eijk (from eq 3), which ranged from -6 to +6 kJ/mol. Many-body energies are often decomposed into terms for polarization, dispersion, induction, and exchange. Such a decomposition can provide physical insight. In our case, however, especially since it has been shown that four- and higher-molecule terms are quite unimportant for water,46–48 it is not useful to perform this kind of decomposition. The threemolecule energies calculated above arise primarily from the cooperative nature of hydrogen bonding, meaning that the strength of a hydrogen bond between two molecules depends on the position of the third molecule. For a trimer configuration where the bridging molecule is both a hydrogen bond donor and acceptor (called type B in Figure 2), the three-molecule energy is generally negative (“cooperative”), while if the bridging molecule is a double donor (type A) or double acceptor (type C), the three-molecule energy is generally positive (“anticooperative”).48,49,122–124 Hydrogen bonding interactions depend to a large extent on intermolecular H-O distances. We

Figure 2. The three types of three-molecule interactions.

Figure 3. Definition of trimer distances. On each molecule a denotes an O atom while b and c label the two H atoms. The three molecules are referred to as 1, 2, and 3.

found a good fit to the ab initio three-molecule energies using a form that involves three exponentials, for interactions of types A, B, and C, respectively:

f(r1, r2) ) A1e-A2(r1+r2) g(r1, r2) ) B1e-B2(r1+r2) h(r1, r2) ) C1e-C2(r1+r2)

(5)

which depend on the relevant H-O distances as shown in Figure 2. For a given trimer composed of three molecules labeled 1, 2, and 3 as in Figure 3, there are 12 relevant H-O distances, and we write the three-molecule energy as a sum of three terms, involving interactions of types A, B, and C, respectively:

∆E123 ) ∆EA123 + ∆EB123 + ∆EC123

(6)

∆EA123 is composed of the 6 terms,

∆EA123 ) f(rb1a2, rc1a3) + f(rb1a3, rc1a2) + f(rb2a1, rc2a3) + f(rb2a3, rc2a1) + f(rb3a1, rc3a2) + f(rb3a2, rc3a1), (7) ∆EB123 is composed of the 24 terms,

8314 J. Phys. Chem. B, Vol. 112, No. 28, 2008

Kumar and Skinner

∆EB123 ) g(rb1a2, rb3a1) + g(rb1a2, rc3a1) + g(rc1a2, rb3a1) + g(rc1a2, rc3a1) + g(rb1a3, rb2a1) + g(rb1a3, rc2a1) + g(rc1a3, rb2a1) + g(rc1a3, rc2a1) + g(rb2a1, rb3a2) + g(rb2a1, rc3a2) + g(rc2a1, rb3a2) + g(rc2a1, rc3a2) + g(rb2a3, rb1a2) + g(rb2a3, rc1a2) + g(rc2a3, rb1a2) + g(rc2a3, rc1a2) + g(rb3a1, rb2a3) + g(rb3a1, rc2a3) + g(rc3a1, rb2a3) + g(rc3a1, rc2a3) + g(rb3a2, rb1a3) + g(rb3a2, rc1a3) + g(rc3a2, rb1a3) + g(rc3a2, rc1a3), (8) and

∆EC123

is composed of the 12 terms:

∆EC123 ) h(rb2a1, rb3a1) + h(rb2a1, rc3a1) + h(rc2a1, rb3a1) +

Figure 4. Scatter plot of the three-molecule fit energies ∆E123 versus a ab initio energies ∆E123 . The diagonal line is a guide to the eye. The root-mean-square deviation is 0.69 kJ/mol.

h(rc2a1, rc3a1) + h(rb1a2, rb3a2) + h(rb1a2, rc3a2) + h(rc1a2, rb3a2) + h(rc1a2, rc3a2) + h(rb1a3, rb2a3) + h(rb1a3, rc2a3) + h(rc1a3, rb2a3) + h(rc1a3, rc2a3), (9) where the distances are between H and O atoms as labeled in Figure 3. To obtain the parameters A1, A2, B1, B2, C1, and C2, we performed a global fit to the ab initio three-molecule energies using a genetic algorithm.125,126 The values for the parameters are listed in Table 1; note that in agreement with the above discussion, A1, C1 > 0, and B1 < 0. Figure 4 shows a scatter plot of the fit three-molecule energies ∆E123 versus those from the ab initio calculations ∆Ea123. The rms deviation of our fit from the ab initio results is 0.69 kJ/mol. C. Simulation of the Liquid. These results for the two- and three-molecule energies would presumably lead to a reasonable approximation to the true Born-Oppenheimer surface for the liquid through eq 2. However, as discussed above, this is not the appropriate potential for a classical simulation of the liquid. Therefore, the idea is to develop an effective classical potential whose form is based on these ab initio results, but whose magnitude is determined empirically. Thus we write:

Eij ) Eijr + λ1∆Eij

(10)

A C B + ∆Eijk ) + λ3∆Eijk ∆Eijk ) λ2(∆Eijk

(11)

λ1, λ2, and λ3 are parameters to be determined by comparing classical simulation results for this potential with experiment. As indicated, for simplicity we decided to use the same parameter λ2 for both types of anticooperative interactions. Note that λ1 ) λ2 ) λ3 ) 1 gives the (approximate) ab initio potential, while λ1 ) λ2 ) λ3 ) 0 recovers the reference (pairwise) potential. In performing an MD simulation with this new potential it is especially efficient if the three-molecule energies can be truncated at some appropriate point. We observe that the threemolecule terms are nearly zero when H-O distances exceed 5 TABLE 1: The parameters used in the three-molecule energies parameter

value

A1 A2 B1 B2 C1 C2

4699.6 kJ/mol 1.88 Å -1 -2152.9 kJ/mol 1.71 Å -1 1312.7 kJ/mol 1.56 Å -1

Å. To effect a smooth cutoff, we multiply each of the exponentials in eq 5 by the product of switching functions s(r1)s(r2), where127

s(r) )

(rf - r)2(rf + 2r - 3rs) (rf - rs)3

, rs e r e rf

(12)

and s(r) ) 1 for r < rs and 0 for r > rf. rs is set to 5.0 Å and rf to 5.2 Å. We checked to see that this modification did not appreciably affect either the fit to the ab initio energies or the total potential energy in the simulations. For several sets of parameters λi, classical molecular dynamics simulations were performed in the NVE ensemble. Our simulation box contained 256 molecules and the size of the box was chosen to give the experimental number density of water at 298 K (3.32 × 1028 m -3).128 We used periodic boundary conditions, and approximations to the Ewald sum were employed to take into account long-range contributions to the electrostatic interactions.129 Rotations were treated using quaternions,130 and the equations of motion were integrated using the leapfrog algorithm131 with a time step of 1 fs. The system was equilibrated to 298 K by rescaling the velocities of the particles. To determine approximate error bars for the various dynamical quantities (see below), we divided the trajectories into five blocks and calculated the standard deviation of the mean. D. Empirical Determination of λi. With three adjustable parameters, we can try to fit three (or more) experimental quantities. We decided to focus on the position of the first peak in the O-O radial distribution function, the coordination number (as defined by integrating the O-O radial distribution function out to 3.36 Å, the position of the first minimum in the experimental function), and the diffusion constant; the experimental values of these quantities are, respectively, 2.75 Å,132 4.47,132 and 2.30 × 10-5 cm 2 /s.133 This fitting is a somewhat cumbersome procedure, since for each set of λi one needs to run a reasonably long simulation to calculate these three quantities. For this reason, a three-dimensional “grid’ approach is untenable. Thus, we developed an alternative iterative approach as follows. Labeling these three quantities as Pj, j ) 1-3, we can consider each of the Pj to be a function of the three λi. Thus

dPj )

∂P

∑ ∂λij dλi

(13)

i 0 0 0 Starting at an initial point (λ1, λ2, λ3), a total of seven simulations at this point and at the points (λ01 ( δλ, λ02, λ03), (λ01, λ02 ( δλ, λ03)

Water Simulation Model

J. Phys. Chem. B, Vol. 112, No. 28, 2008 8315

TABLE 2: The Position (R) of the First Maximum in the O-O Radial Distribution Function, the Coordination Number (C), and the Diffusion Constant (D) for the SPC/E Model, the TIP4P Model, and our New Model together with the Experimental Values model

R (Å)

C

SPC/E TIP4P new model experiment

2.76 2.76 2.76 2.75132

4.55 4.56 4.47 4.47132

D (10-5 cm2/s) 2.49 ( 0.04 3.44 ( 0.07 2.36 ( 0.05 2.30133

and (λ01, λ02, λ03 ( δλ) are run. At each point, the quantities Pj are calculated, leading to estimates for ∂Pj/∂λi. Next, we consider a finite difference approximation to the above:

∆Pj )

∂P

∑ ∂λij ∆λi

(14)

i

Taking ∆Pj to be the difference between the experimental value of Pj and the calculated value at the point λ0i , this equation is then solved for ∆λi. This leads to new trial values λ1i ) λ0i + ∆λi, and the process is iterated until convergence is obtained. To implement this procedure in principle one can start at any value of λi; reasonable choices are λi ) 0 (the TIP4P reference potential) or λi ) 1 (the ab initio potential). We chose the latter and took δλ ) 0.1. This procedure led to convergence in only two iterations to λ1 ) 0.38, λ2 ) 1.07, and λ3 ) 0.92. As a result of this procedure, the simulation results with these parameters lead to values for these three quantities that are virtually identical to experiment. These values, together with those for the TIP4P and SPC/E models, are shown in Table 2. Especially, for the diffusion constant, the new model provides a considerable improvement over the TIP4P model.

Figure 5. O-O radial distribution function for liquid water at 298 K for the SPC/E model, the TIP4P model, the new model, and the experiment.132

III. Results In the previous section we used three specific experimental quantities for liquid water to parametrize our model. In this section, we calculate several other quantities for this new model, in each case comparing with experiment. First, we find that the internal energy is -40.7 kJ/mol, reasonably close to the experimental value of -41.5 kJ/mol. These can be compared with -41.4 kJ/mol, the value for both the SPC/E and TIP4P models. Second, we can compare the full O-O radial distribution function (not just the distance at the maximum and the coordination number). Our results for the new model are compared to those from experiment132 and the TIP4P and SPC/E models in Figure 5. One sees that the three models yield very similar results, which all are slightly more structured (higher first peak) than experiment. We can make similar comparisons for the H-H and O-H radial distribution functions, as shown in Figures 6 and 7. In both cases, the results of the three simulation models are similar; they agree quite well with experiment for gHH(r), while for gOH(r), the simulation models all produce a first (intermolecular) peak that is too high compared with experiment. Third, we can calculate certain rotational correlation times. The relevant rotational time-correlation function for moleculefixed unit vector uˆ is given by 〈P2(uˆ(0) · uˆ(t))〉 (P2 is the second Legendre polynomial), and the rotational correlation time is the integral from 0 to ∞ of this time-correlation function. To obtain this correlation time for the different models, the rotational timecorrelation functions were calculated to 10 ps and then extrapolated to longer times using an exponential fit. It is of

Figure 6. H-H radial distribution function for liquid water at 298 K for the SPC/E model, the TIP4P model, the new model, and the experiment.132

Figure 7. O-H radial distribution function for liquid water at 298 K for the SPC/E model, the TIP4P model, the new model, and the experiment.132

particular interest to calculate the rotational correlation time of the molecule-fixed intramolecular H-H unit vector, as this is known to be substantially too fast for the TIP4P model (compared to experiment). The values, τHH, for the three models are given in Table 3. One sees that the value for the new model (2.11 ps) is 50% longer than that for the TIP4P model (1.41

8316 J. Phys. Chem. B, Vol. 112, No. 28, 2008

Kumar and Skinner

TABLE 3: Rotational Correlation Times τHH, τOD, and τp, Corresponding to the Different Molecule-Fixed Unit Vectors in H2O, D2O, and D17 2 O, respectively, for Different Simulation Models at 298 K model SPC/E TIP4P New model Experiment

τHH (ps) 1.92 ( 0.03 1.41 ( 0.02 2.11 ( 0.05 2.36134

τOD (ps)

τp (ps)

1.95 ( 0.03 1.43 ( 0.03 2.10 ( 0.06 2.4 135

1.28 ( 0.03 1.04 ( 0.02 1.48 ( 0.02 1.8 135

ps) and is in substantially better agreement with the experiment134 (2.36 ps). Note that this latter value has been obtained by interpolation from the experimental results at 283 and 303 K using an Arrhenius form. Rotational correlation times for other molecule-fixed unit vectors have also been measured.135 Specifically, deuterium quadrupole relaxation experiments on D2O give a value for the correlation time of the OD bond vector of τOD ) 2.4 ps, and 17O quadrupolar relaxation experiments on dilute D17O in D16O 2 2 give a value for the correlation time for the out-of-plane (perpendicular) unit vector of τp ) 1.8 ps. To make a reasonable comparison with these experiments, we performed simulations of neat D2O and neat D17 2 O, but at the experimental H2O density, and calculated τOD and τp for these two systems, respectively. The theoretical values for the SPC/E, TIP4P, and new models are given in Table 3. Again we see that the values for the new model are somewhat longer than for both SPC/E and TIP4P models (especially the latter) and are also are in better agreement with experiment. In addition, the significant difference in the OD and out-of-plane correlation times provides further theoretical support for the idea that water rotation is anisotropic.135 IV. Concluding Remarks We have developed a new classical water simulation model that has explicit two- and three-molecule interactions. The model is formulated as a perturbation to the well-known two-molecule interaction TIP4P model. The forms of the correction to the TIP4P two-molecule term and explicit three-molecule terms have been determined by ab initio calculations on dimers and trimers selected from liquid-state configurations. The resulting potential surface presumably is a reasonable approximation to the true Born-Oppenheimer potential, in which case in order to describe the properties of liquid water one would have to perform nuclear quantum dynamics simulations. Instead, we are interested in developing a model for use with classical simulations, and so we determined the strength of these perturbations empirically by comparing classical simulation results to experiment. The resulting classical model is in better agreement with experiment than the TIP4P (and also the SPC/E) model for important dynamical quantities such as the diffusion constant and various rotational correlation times. Because of the shortrange nature of the three-molecule interactions, the model is not much slower to simulate than a simple two-molecule model (our code runs 1.7 times slower than that for TIP4P, but that is without making any attempt to optimize it). As a further test of the model it would be useful and interesting to calculate IR and Raman spectra and also ultrafast vibrational spectroscopy observables like echo peak shifts and 2DIR spectra, using theoretical techniques that we have developed14,136,137 and compare with the experiment.138,139 It would also be interesting to perform an analysis of hydrogenbonding statics and dynamics, similar to what we have recently done for SPC/E water,140 to shed further light on the important issue of the hydrogen-bonding network in liquid water.141

Our hope is that this new simulation model, by way of its explicit three-molecule interactions, will be applicable to a wide variety of physical situations, from liquid water at higher and lower temperatures to supercritical water to the various phases of ice. Thus our plan will be to calculate the water phase diagram with this model,15–18 in order to assess this broader applicability. In addition, we will calculate the second and higher virial coefficients and compare with experiment (the simple twomolecule interaction classical models do not do well in this regard).142 If these calculations support the broader validity of the model, then it can be used in a variety of applications, especially in systems involving substantial heterogeneities, like the liquid/vapor interface,19–30 ionic and other aqueous solutions,31,32 and reverse micelles.33–37 Acknowledgment. The authors thank the National Science Foundation for support of this work through Grants CHE0446666 and CHE-0750307. J.L.S. thanks Frank Weinhold for helpful discussions about quantum chemistry in general and about cooperativity and hydrogen bonding. References and Notes (1) Guillot, B. J. Mol. Liq. 2002, 101, 219. (2) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B.,. Ed.; Reidel: Dordrecht, 1981. (3) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (4) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (5) Reimers, J. R.; Watts, R. O.; Klein, M. L. Chem. Phys. 1982, 64, 95. (6) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665. (7) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (8) Rick, S. W. J. Chem. Phys. 2004, 120, 6085. (9) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (10) Abascal, J. L. F.; Sanz, E.; Ferna`ndez, R. G.; Vega, C. J. Chem. Phys. 2005, 122, 234511. (11) Amira, S.; Spångberg, D.; Hermansson, K. Chem. Phys. 2004, 303, 327. (12) Wu, Y.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (13) Lawrence, C. P.; Skinner, J. L. Chem. Phys. Lett. 2003, 372, 842. (14) Schmidt, J. R.; Roberts, S. T.; Loparo, J. J.; Tokmakoff, A.; Fayer, M. D.; Skinner, J. L. Chem. Phys. 2007, 341, 143. (15) Vega, C.; Abascal, J. L. F.; Sanz, E.; MacDowell, L. G.; McBride, C. J. Phys.: Condens. Matter 2005, 17, 3283. (16) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. (17) Vega, C.; Abascal, J. L. F.; Nezbeda, I. J. Chem. Phys. 2006, 125, 034503. (18) Baranyai, A.; Barto´k, A.; Chialvo, A. A. J. Chem. Phys. 2006, 124, 074507. (19) Shen, Y. R.; Ostroverkhov, V. Chem. ReV. 2006, 106, 1140. (20) Brown, M. G.; Raymond, E. A.; Allen, H. C.; Scatena, L. F.; Richmond, G. L. J. Phys. Chem. A 2000, 104, 10220. (21) Raymond, E. A.; Tarbuck, T. L.; Brown, M. G.; Richmond, G. L. J. Phys. Chem. B 2003, 107, 546. (22) Shultz, M. J.; Baldelli, S.; Schnitzer, C.; Simonelli, D. J. Phys. Chem. B 2002, 106, 5313. (23) Liu, D.; Ma, G.; Levering, L. M.; Allen, H. C. J. Phys. Chem. B 2004, 108, 2252. (24) Buch, V. J. Phys. Chem. B 2005, 109, 17771. (25) Benjamin, I. Phys. ReV. Lett. 1994, 73, 2083. (26) Morita, A.; Hynes, J. T. Chem. Phys. 2000, 258, 371. (27) Morita, A.; Hynes, J. T. J. Phys. Chem. B 2002, 106, 673. (28) Perry, A.; Neipert, C.; Space, B.; Moore, P. B. Chem. ReV. 2006, 106, 1234. (29) Brown, E.; Mucha, M.; Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2005, 109, 7934. (30) Walker, D. S.; Hore, D. K.; Richmond, G. L. J. Phys. Chem. B 2006, 110, 20451. (31) Park, S.; Fayer, M. D. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 16731.

Water Simulation Model (32) Smith, J. D.; Saykally, R. J.; Geissler, P. L. J. Am. Chem. Soc. 2007, 129, 13847. (33) Faeder, J.; Ladanyi, B. M. J. Phys. Chem. B 2000, 104, 1033. (34) Harpham, M. R.; Ladanyi, B. M.; Levinger, N. E.; Herwig, K. W. J. Chem. Phys. 2004, 121, 7855. (35) Harpham, M. R.; Ladanyi, B. M.; Levinger, N. E. J. Phys. Chem. B 2005, 109, 16891. (36) Piletic, I. R.; Moilanen, D. E.; Spry, D. B.; Levinger, N. E.; Fayer, M. D. J. Phys. Chem. A 2006, 110, 4985. (37) Spry, D. B.; Goun, A.; Glusac, K.; Moilanen, D. E.; Fayer, M. D. J. Am. Chem. Soc. 2007, 129, 8122. (38) Kuharski, R. A.; Rossky, P. J. J. Chem. Phys. 1985, 82, 5164. (39) Poulsen, J. A.; Nyman, G.; Rossky, P. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6709. (40) Poulsen, J. A.; Nyman, G.; Rossky, P. J. J. Chem. Theory Comput. 2006, 2, 1482. (41) Stern, H. A.; Berne, B. J. J. Chem. Phys. 2001, 115, 7622. (42) Fanourgakis, G. S.; Schenter, G. K.; Xantheas, S. S. J. Chem. Phys. 2006, 125, 141102. (43) Paesani, F.; Zhang, W.; Cheatum, T. E.; Voth, G. A. J. Chem. Phys. 2006, 125, 184507. (44) Paesani, F.; Iuchi, S.; Voth, G. A. J. Chem. Phys. 2007, 127, 074506. (45) Tsai, C. J.; Jordan, K. D. J. Phys. Chem. 1993, 97, 5208. (46) Pedulla, J. M.; Vila, F.; Jordan, K. D. J. Chem. Phys. 1996, 105, 11091. (47) Hodges, M. P.; Stone, A. J.; Xantheas, S. S. J. Phys. Chem. A 1997, 101, 9163. (48) Ojama¨e, L.; Hermansson, K. J. Phys. Chem. 1994, 98, 4271. (49) Xantheas, S. S. Chem. Phys. 2000, 258, 225. (50) Caldwell, J.; Dang, L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 9144. (51) Dang, L. X. J. Chem. Phys. 1992, 97, 2659. (52) Dang, L. X.; Chang, T.-M. J. Chem. Phys. 1997, 106, 8149. (53) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (54) Lamoureux, G.; MacKerell, A. D.; Roux, B. J. Chem. Phys. 2003, 119, 5185. (55) Chen, B.; Xing, J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2391. (56) Baranyai, A.; Bartok, A. J. Chem. Phys. 2007, 126, 184508. (57) Stern, H. A.; Rittner, F.; Berne, B. J.; Friesner, R. A. J. Chem. Phys. 2001, 115, 2237. (58) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933. (59) Paricaud, P.; Predota, M.; Chialvo, A. A.; Cummings, P. T. J. Chem. Phys. 2005, 122, 244511. (60) Iuchi, S.; Izvekov, S.; Voth, G. A. J. Chem. Phys. 2007, 126, 124505. (61) Bursulaya, B. D.; Kim, H. J. J. Chem. Phys. 1998, 108, 3277. (62) Bursulaya, B. D.; Jeon, J.; Zichi, D. J.; Kim, H. J. J. Chem. Phys. 1998, 108, 3286. (63) Lefohn, A. E.; Ovchinnikov, M.; Voth, G. A. J. Phys. Chem. B 2001, 105, 6628. (64) Jeon, J.; Lefohn, A. E.; Voth, G. A. J. Chem. Phys. 2003, 118, 7504. (65) Kuwajima, S.; Warshel, A. J. Phys. Chem. 1990, 94, 460. (66) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (67) Piquemal, J.-P.; Chelli, R.; Procacci, P.; Gresh, N. J. Phys. Chem. A 2007, 111, 8170. (68) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 5, 2471. (69) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 105, 1142. (70) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2002, 116, 10372. (71) Chen, B.; Ivanov, I.; Klein, M. L.; Parrinello, M. Phys. ReV. Lett. 2003, 91, 215503. (72) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 120, 300. (73) Kuo, I.-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I. J. Chem. Theory Comput. 2006, 2, 1274. (74) Lee, H.-S.; Tuckerman, M. E. J. Chem. Phys. 2006, 125, 154507. (75) Lee, H.-S.; Tuckerman, M. E. J. Chem. Phys. 2007, 126, 164501. (76) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 3572. (77) Schwegler, E.; Grossman, J. C.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 121, 5400. (78) Mantz, Y. A.; Chen, B.; Martyna, G. J. J. Phys. Chem. B 2006, 110, 3540. (79) Kuo, I.-F. W.; Mundy, C. J. Science 2004, 303, 658. (80) Todorova, T.; Seitsonen, A. P.; Hutter, J.; Kuo, I.-F. W.; Mundy, C. J. J. Phys. Chem. B 2006, 110, 3685. (81) Schwegler, E.; Galli, G.; Gygi, F. Phys. ReV. Lett. 2000, 84, 2429. (82) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515.

J. Phys. Chem. B, Vol. 112, No. 28, 2008 8317 (83) Matsuoka, O.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976, 64, 1351. (84) Lie, G. C.; Clementi, E.; Yoshimine, M. J. Chem. Phys. 1976, 64, 2314. (85) Wojcik, M.; Clementi, E. J. Chem. Phys. 1986, 84, 5970. (86) Niesar, U.; Corongiu, G.; Clementi, E.; Kneller, G. R.; Bhattacharya, D. K. J. Phys. Chem. 1990, 94, 7949. (87) Burnham, C. J.; Li, J. C.; Xantheas, S. S. J. Chem. Phys. 1999, 110, 4566. (88) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115. (89) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1479. (90) Millot, C.; Stone, A. J. Mol. Phys. 1992, 77, 439. (91) Saykally, R. J.; Blake, G. A. Science 1993, 259, 1570. (92) Fellers, R. S.; Leforestier, C.; Braly, L. B.; Brown, M. G.; Saykally, R. J. Science 1999, 284, 945. (93) Goldman, N.; Leforestier, C.; Saykally, R. J. Philos. Trans. R. Soc. London, Ser. A 2004, 363, 493. (94) Mas, E. M.; Szalewicz, K.; Bukowski, R.; Jeziorski, B. J. Chem. Phys. 1997, 107, 4207. (95) Mas, E. M.; Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; Wormer, P. E. S.; van der Avoird, A. J. Chem. Phys. 2000, 113, 6687. (96) Mas, E. M.; Bukowski, R.; Szalewicz, K. J. Chem. Phys. 2003, 118, 4404. (97) Mas, E. M.; Bukowski, R.; Szalewicz, K. J. Chem. Phys. 2003, 118, 4386. (98) Bukowski, R.; Szalewicz, K.; Groenenboom, G.; van der Avoird, A. J. Chem. Phys. 2006, 125, 044301. (99) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Science 2007, 315, 1249. (100) Yoon, B. J.; Morokuma, K.; Davidson, E. R. J. Chem. Phys. 1985, 83, 1223. (101) Patkowski, K.; Szalewicz, K.; Jeziorski, B. J. Chem. Phys. 2006, 125, 154107. (102) Groenenboom, G. C.; Mas, E. M.; Bukowski, R.; Szalewicz, K.; Wormer, P. E. S.; van der Avoird, A. Phys. ReV. Lett. 2000, 84, 4072. (103) Keutsch, F. N.; Saykally, R. J. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 10533. (104) Goldman, N.; Fellers, R. S.; Brown, M. G.; Braly, L. B.; Keoshian, C. J.; Leforestier, C.; Saykally, R. J. J. Chem. Phys. 2002, 116, 10148. (105) Leforestier, C.; Gatti, F.; Fellers, R. S.; Saykally, R. J. J. Chem. Phys. 2002, 117, 8710. (106) Harker, H. A.; Keutsch, F. N.; Leforestier, C.; Scribano, Y.; Han, J.-X.; Saykally, R. J. Mol. Phys. 2007, 105, 497. (107) Harker, H. A.; Keutsch, F. N.; Leforestier, C.; Scribano, Y.; Han, J.-X.; Saykally, R. J. Mol. Phys. 2007, 105, 513. (108) Fanourgakis, G. S.; Xantheas, S. S. J. Phys. Chem. A 2006, 110, 4100. (109) Netzloff, H. M.; Gordon, M. S. J. Chem. Phys. 2004, 121, 2711. (110) Wallqvist, A.; Karlstro¨m, G. Chem. Scr. 1989, 29, 131. (111) Wallqvist, A.; Ahlstro¨m, P.; Karlstro¨m, G. J. Phys. Chem. 1990, 94, 1649. (112) Engkvist, O.; Forsberg, N.; Schu¨tz, M.; Karlstro¨m, G. Mol. Phys. 1997, 90, 277. (113) Spångberg, D.; Hermansson, K. J. Chem. Phys. 2003, 119, 7263. (114) Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. J. Chem. Phys. 2004, 120, 10896. (115) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2005, 123, 134105. (116) Gaussian 03. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al., G. E. S. Gaussian Inc.: Pittsburgh PA, 2003. (117) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (118) Liu, B.; McLean, A. D. J. Chem. Phys. 1973, 59, 4557. (119) Xantheas, S. S. J. Chem. Phys. 1996, 104, 8821. (120) Xantheas, S. S.; Burnham, C. J.; Harrison, R. J. J. Chem. Phys. 2002, 116, 1493. (121) Halkier, A.; Koch, H.; Jørgensen, P.; Christiansen, O.; Nielsen, I. M. B.; Helgaker, T. Theor. Chem. Acc. 1997, 97, 150. (122) Reed, A. E.; Weinhold, F. J. Chem. Phys. 1983, 78, 4066. (123) Reed, A. E.; Weinhold, F.; Curtiss, L. A.; Pochatko, D. J. J. Chem. Phys. 1986, 84, 5687. (124) Ohno, K.; Okimura, M.; Akai, N.; Katsumoto, Y. Phys. Chem. Chem. Phys. 2005, 7, 3005. (125) Michalewicz, Z. Genetic Algorithms + Data Structures ) EVolution Programs, 3rd ed; Springer-Verlag: Berlin, 1996. (126) Goldberg, D. E. Genetic Algorithms in Search, Optimization, and Machine Learning; Addison-Wesley: Reading, MA, 1989. (127) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (128) Franks, F. The physics and physical chemistry of water; Plenum Press: New York, 1972; Vol. 1. (129) Adams, D. J.; Dubey, G. S. J. Comput. Phys. 1987, 72, 156. (130) Svanberg, M. Mol. Phys. 1997, 92, 1085. (131) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987.

8318 J. Phys. Chem. B, Vol. 112, No. 28, 2008 (132) Soper, A. K. Chem. Phys. 2000, 258, 121. (133) Krynicki, K.; Green, C. D.; Sawyer, D. W. Faraday Discuss. Chem. Soc. 1978, 66, 199. (134) Jonas, J.; DeFries, T.; Wilbur, D. J. J. Chem. Phys. 1976, 65, 582. (135) Ropp, J.; Lawrence, C.; Farrar, T. C.; Skinner, J. L. J. Am. Chem. Soc. 2001, 123, 8047. (136) Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2004, 120, 8107. (137) Auer, B.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14215. (138) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698.

Kumar and Skinner (139) Asbury, J. B.; Steinel, T.; Stromberg, C.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. J. Phys. Chem. A 2004, 108, 1107. (140) Kumar, R.; Schmidt, J. R.; Skinner, J. L. J. Chem. Phys. 2007, 126, 204107. (141) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Na¨slund, L. Å.; Hirsch, T. K.; Ojama¨e, L.; Glatzel, P.; Pettersson, L. G. M.; Nilsson, A. Science 2004, 304, 995. (142) Benjamin, K. M.; Singh, J. K.; Schultz, A. J.; Kofke, D. A. J. Phys. Chem. B 2007, 111, 11463.

JP8009468