The Properties of Water: Insights from Quantum Simulations - The

5702

J. Phys. Chem. B 2009, 113, 5702–5719

CENTENNIAL FEATURE ARTICLE The Properties of Water: Insights from Quantum Simulations† Francesco Paesani and Gregory A. Voth* Center for Biophysical Modeling and Simulation and Department of Chemistry, UniVersity of Utah, 315 South 1400 East Room 2020, Salt Lake City, Utah 84112-0850 ReceiVed: December 2, 2008

The properties of water play a central role in many phenomena of relevance to different areas of science, including physics, chemistry, biology, geology, and climate research. Although well studied for decades, the behavior of water under different conditions and in different environments still remains mysterious and often surprising. In this article, various efforts aimed at providing a comprehensive representation of the water properties at a molecular level through computer modeling and simulation will be described. In particular, the unique role played by the hydrogen-bond network will be examined, first in liquid water, then in the solvation of model biological compounds, and finally in ice, especially highlighting the important effects related to the quantization of the nuclear motion. 1. Introduction Given the fundamental role played by water in a multitude of phenomena, it is not surprising that many theoretical and experimental studies have been devoted in the past to the development of a detailed understanding at a molecular level of its unique properties. From a theoretical point of view, numerous water models have been developed, ranging from very simple ones based on atomic partial charges and rigid bonds1-9 to other more sophisticated ones that include molecular flexibility10-16 and explicit electronic polarization effects.17-34 When employed in computer simulations, each of these models has shown a certain degree of success in reproducing at least some of the properties of liquid water and ice. The theoretical investigations have been accompanied by many experiments using X-ray and neutron scattering, NMR, Raman and infrared (IR) spectroscopy, etc.35 Combined together, these studies certainly make water the molecular species for which the greatest number of data has been collected. Unfortunately, such efforts have not always been followed by a corresponding improvement in our understanding of water properties. As a proof of this, only a few years ago, the exact determination of the structure of water was listed among the most outstanding problems in science.36 In this context, a tenuous consensus on the average structure of liquid water has been achieved only in the past few years from independent X-ray37 and neutron scattering38 experiments, while NMR measurements have been used to determine the angular distribution of the hydrogen bonds as a function of the temperature.39 However, the generally accepted picture of 4-fold † 2008 marked the Centennial of the American Chemical Society’s Division of Physical Chemistry. To celebrate and to highlight the field of physical chemistry from both historical and future perspectives, The Journal of Physical Chemistry is publishing a special series of Centennial Feature Articles. These articles are invited contributions from current and former officers and members of the Physical Chemistry Division Executive Committee and from J. Phys. Chem. Senior Editors. * Author to whom correspondence should be addressed. E-mail: [email protected].

coordinated water molecules in a tetrahedral arrangement has also been challenged. From an analysis of X-ray absorption spectra measured for liquid water and ice, complemented with ab initio calculations for an ensemble of representative water clusters, it has been suggested that the coordination number in the liquid is actually closer to two with one donor and one acceptor hydrogen bond for each molecule.40 This proposed structural organization with hydrogen-bonded chains or large rings of water molecules embedded in a weakly hydrogenbonded disordered network is still highly controversial.41-49 On the other hand, infrared spectra of liquid water measured at different temperatures have demonstrated the importance of librations (or hindered rotations) for the formation of a dynamic hydrogen-bond (H-bond) network,35 while the presence of isosbestic points in the Raman spectra suggested the existence of two distinct types of H2O environments.50 This hypothesis found support from earlier ultrafast spectroscopic studies that suggested the existence of H2O molecules having two different relaxation rates.51 However, recent measurements of the temperature dependence of Raman spectra combined with molecular dynamics (MD) simulations have shown that the experimental data can also be interpreted in terms of a continuous distribution of structures.52 Importantly, the enthalpy differences related to the possible water structures have been estimated from experiment to be approximately equal to 2.6 kcal mol-1 (900 cm-1) from Raman spectroscopy,53 or to 1.5 kcal mol-1 (500 cm-1) from X-ray absorption spectroscopy.50 These values, which lie in the range of the energies required for the librational motion, strongly suggest that water molecules in both fundamental and first excited librational states coexist in the liquid where they display large rotational amplitudes. As a consequence, nuclear quantum effects are expected to play a nontrivial role in determining the dynamics of the H-bond network. In this regard, recent narrowband two-color pump-probe measurements of the orientational dynamics of dilute HOD molecules in D2O showed a different frequency behavior from that observed for HOD molecules in H2O.54,55 In order to explain this difference, it was

10.1021/jp810590c CCC: $40.75  2009 American Chemical Society Published on Web 03/19/2009

Centennial Feature Article

Francesco Paesani received his Ph.D. in Theoretical Physical Chemistry from the University of Rome “La Sapienza”, Italy. His present research includes the development of theoretical and computational methodologies for the investigation of structural, thermodynamic, and dynamical properties of many-body systems at a quantum-mechanical level, including multiphase complex systems such as atmospherically relevant interfaces. Dr. Paesani is currently a postdoctoral researcher in the group of Professor Gregory Voth at the University of Utah. He was previously a postdoctoral researcher in the group of Birgitta Whaley at the University of California, Berkeley.

Gregory A. Voth is Distinguished Professor of Chemistry and Director of the Center for Biophysical Modeling and Simulation at the University of Utah. He has authored or co-authored approximately 300 scientific publications and mentored over 130 postdoctoral fellows and graduate students. His research interests include the theory and simulation of charge transport, solvation, and multiscale phenomena in liquid state, biomolecular, and materials systems. He was Chair of the American Chemical Society Division of Physical Chemistry in 2008.

proposed that tunneling could selectively speed up the reorientation of the most weakly hydrogen-bonded OH groups.55 By contrast, no frequency dependence was found in similar IR pump-probe experiments for the HOD:D2O system with 45 fs pulses, which had sufficient bandwidth to span the entire OH absorption band.56 More recent studies of the HOD:H2O system with the pump tuned to different frequencies of the absorption band also found a frequency dependence in the OD orientational dynamics for delays 1 thermostats to each degree of freedom of the system.98 Unfortunately, a direct solution of eq 6 using standard MD algorithms is doomed to failure, as established in ref 99. In fact, as P becomes large, the effective force constant between the beads grows linearly with P, while the contribution from V is attenuated by a factor of 1/P. Consequently, the harmonic term dominates, which causes the MD trajectories to not explore the


J. Phys. Chem. B, Vol. 113, No. 17, 2009 5705

full available phase space. In order to solve this problem, it is necessary to uncouple the harmonic term in the Hamiltonian of eq 7, which can be accomplished by defining a more convenient coordinate system. In our quantum simulations, we adopted the so-called normal mode representation of the path-integrals.93 Within this scheme, the Cartesian coordinates of the ring-polymer beads are transformed in normal modes through the following unitary transformation P

qn )

∑

1 Unixj √P j)1

(9)

P

∑ Ujnqn

xj ) √P

(10)

The centroid molecular dynamics formalism draws upon the prescription of distribution functions, in which the exact quantum expressions are cast into a phase space representation leading to a classical-like interpretation of the variables of interest. Within CMD, the real-time dynamics of a quantum particle is described by the following classical-like equation of motion

mq¨c(t) ) 〈F(qc)〉c

where qc is the centroid position of the quantum particle, while the right-hand side corresponds to the centroid force. The latter is defined by the following constrained path-integral average

〈F(qc)〉c ) P

[

i)1

where

1 2πjn 2πjn Ujn ) cos - sin P P √P

(

∫

)

)

(

]

)

P

∑

]

(15)

In eq 15, q˜ is the centroid variable corresponding to the center of mass of the ring-polymer in a path-integral representation of the quantum particle, i.e., P

q˜ ) 2πn 1 - cos P

i)1

1 dx1...dxP δ(qc - q˜) exp -β V(xi) P i)1

(11)

Clearly, the unitary matrix U diagonalizes the harmonic term in eq 7, with the eigenvalues being the squared angular frequencies of the nth normal mode

2mωP2

[

P

∫ dx1...dxP δ(qc - q˜)[ ∑ fi] exp -β ∑ P1 V(xi)

n)1

ωn2

(14)

(12)

∑

1 x P i)1 i

(16)

and fi is the force acting on each bead of the ring-polymer Within this representation, the Hamiltonian in eq 7 becomes P

H)

∑

n)1

[

]

P pn2 1 1 + ωn2qn2 + V(√P Ujnqn) 2µ˜ n 2 P n)1

∑

P

∂Φ({xi}) 1 fi ) P i)1 ∂xi

∑

(13)

where µ˜ n ) mωn2 is the fictitious mass of the nth normal mode. After inclusion of an appropriate thermostatting scheme,98 the equations of motion for the normal modes derived from the Hamiltonian in eq 13 can be efficiently propagated using standard MD algorithms resulting in the so-called normal mode path-integral MD (NMPIMD) method. Alternative schemes, such as the staging transformation,100 have also been proposed for decoupling the harmonic term in eq 7. While (numerically) exact quantum thermodynamic and structural properties of water can be computed using the NMPIMD method described above, the situation is much more complicated for quantum dynamical properties. The main reason for this is that an exact description of quantum dynamics in the condensed phase requires the solution of the quantum Liouville equation, which remains one of the most challenging problems in theoretical physical chemistry. Several approximate methods have been developed for describing quantum dynamics in the condensed phase (e.g., centroid molecular dynamics,101-111 semiclassicalapproaches,112-121 quantummodecouplingtheory,122,123 the Feynman-Kleinert linearized path-integral method,124-129 ring-polymer molecular dynamics,130-139 and the effective potential analytic continuation method140,141). In the following, the centroid molecular dynamics (CMD) method that has been extensively developed by our group and employed to study several quantum systems will be briefly reviewed.

(17)

In eq 17, Φ({xi}) is the potential term of the Hamiltonian in eq 7, i.e., P

Φ({xi}) )

∑ [ 21 mωP2(xi - xi+1)2 + P1 V(xi)]

(18)

i)1

Among the different algorithms that have been proposed to propagate the centroid equation of motion in eq 14,104 the adiabatic centroid molecular dynamics (ACMD) scheme can be efficiently implemented within MD. ACMD makes use of the normal mode transformation in eqs 9-11, from which it is possible to show that the (zero-frequency) Pth mode corresponds to the centroid position of the quantum particle, i.e., qP ) qc. In order to take full advantage of the normal mode transformation, the centroid mass is then set equal to that of the quantum particle, µ˜ P ) m, while the mass of each nonzero-frequency mode is set to µ˜ n ) γ2mωn2, where γ is called the adiabaticity parameter.104 The latter has to be chosen small enough in order to ensure a sufficiently large separation in time scales between the centroid and all noncentroid motions. The equations of motion in the ACMD method are then derived from the Hamiltonian in eq 13 as in the case of NMPIMD. However, since the centroid variable must strictly follow a Newtonian dynamics, only the motion of the nonzero-frequency modes must be thermostatted. Within ACMD, the centroid thus moves on

5706 J. Phys. Chem. B, Vol. 113, No. 17, 2009 the potential of mean force generated by the nonzero-frequency modes according to the adiabatic theorem.104 It has been demonstrated that the time propagation of the centroid position and momentum provides an efficient means to compute approximate quantum time correlation functions involving linear operators.106,107 In the case of nonlinear operators, the situation is more complicated, which has led to the development of different numerical schemes.108,109 Recently, a new approach, which combines CMD with the maximum entropy analytic continuation method, has been shown to provide remarkably good agreement with the exact quantum results for several model systems.111 In order to investigate the behavior of water, both the NMPIMD and CMD methods described above must be combined with an accurate representation of the underlying water potential energy function. The first quantum calculations on water were performed in the mid-1980s, using PIMC with a rigid water model to analyze the impact of quantum fluctuations on the liquid structure.142,143 Later, the effects of nuclear quantization on the liquid dynamics were also investigated. CMD was initially applied to study the properties of water at ambient conditions,16 while the less sophisticated Feynman-Hibbs effective potential approach was employed to compute the selfdiffusion coefficient of light (H2O) and heavy (D2O) water over a rather wide range of temperatures.144 The CMD method was also used with a rigid-body water model in several studies of the orientational and translational motions.145-147 More recently, the static and dynamical properties of a flexible water model have been computed using the Feynman-Kleinert linearized path-integral method,127 while the ring-polymer molecular dynamics approach with a rigid-body water model has been applied to the study of quantum diffusion and relaxation processes.148 However, the water models employed in these quantum simulations were originally parametrized using classical mechanics and, as a consequence, the agreement between the calculated and experimental results was found, in general, not to be satisfactory. To solve this problem, empirical quantum force fields were also developed, which showed an improved agreement with the experimental results when employed in PIMD and CMD calculations of the water properties at ambient conditions.15,16 A more reliable description of the underlying interactions is clearly needed in order to achieve a realistic modeling of water with a proper accounting of the nuclear quantum effects. In this context, Car-Parrinello molecular dynamics (CPMD) provides an effective scheme to compute the molecular interactions “on the fly” using density functional theory.149 However, it has been shown that some of the water properties are sensitive to the density functional model chosen to compute the molecular interactions, with different functional models and basis sets providing significantly different results.150-153 Moreover, due to the associated computational cost, quantum CPMD simulations have been limited to the study of structural properties,43,154 and no quantum dynamical calculations have been reported so far. Recently, the proton momentum distribution in liquid water has also been analyzed using an open path-integral approach combined with Car-Parrinello molecular dynamics.155 Polarizable force fields fitted to high level ab initio data offer a promising alternative for the study of the behavior of water.19,25,26,156-160 Such force fields are in principle transferable and, consequently, well suited for a detailed investigation of the water properties under different conditions and in different environments. Furthermore, they are less computationally demanding than the electronic structure calculations required

Paesani and Voth by the CPMD method, and they can thus be used in longer simulations with a larger number of water molecules, which can provide a more statistically accurate description of the liquid phase. Quantum simulations have recently been reported with the polarizable MCDHO,161 QMPFF2,30 TTM2.1-F,162-164 TTM3F,25 and TTM4-F26 water models. 3. Structural Properties of Liquid Water As mentioned in the Introduction, the generally accepted picture of liquid water in which each molecule is, on average, 4-fold coordinated accepting and donating two H-bonds has recently been challenged.40 However, the proposal of a different structural organization consisting of chains and large rings of hydrogen-bonded molecules has not found support from any computer simulations so far.43,44,47,163,165 By contrast, it has been shown that thermodynamic and dielectric properties computed with water models, which were parametrized ad hoc to reproduce a local bonding environment of two H-bonds, are in poor agreement with the available experimental data.166 Although a great deal of information on the structure of liquid water at a molecular level has been obtained from classical MD simulations, X-ray diffraction measurements have shown that the difference in the structure factors between liquid H2O and D2O displays a clear variation between 268 and 318 K,62,63 which suggests that nuclear quantum effects should not be neglected in any truly realistic model of liquid water. In this context, our contributions toward a development of such a model have been via an extensive investigation of its properties through nuclear quantum simulations using NMPIMD and ACMD with ab initiobased, polarizable force fields, such as the TTM2.1-F24 and TTM3-F25 models. It is important to note that, since polarization effects are explicitly included in the parametrization of these models, the atomic charges of each water molecule instantaneously vary in response to any change of the local environment. As a consequence, both models do not assume any a priori structure of the liquid but are instead capable of capturing the occurrence of possible asymmetries in the H-bond network. This, in combination with the absence of implicit nuclear quantum effects in their parametrization, makes such models quite suitable for a comprehensive analysis of the water properties through nuclear quantum simulations. The temperature dependence of the quantum effects on the structure of liquid H2O was first investigated over the interval from 263.15 to 373.15 K using the TTM2.1-F model. For this purpose, the intensity of X-ray scattering was calculated as

I(Q) )

∑ xi xj fi(Q) fj(Q) ij

sin(QRij) + QRij

∑ xi xj fi(Q) fj(Q) Hij(Q)

(19)

i 3 Å, while substantial differences appear with the most recent data of ref 64. The NMPIMD calculations with the TTM2.1-F model also indicate that the intermolecular HH distances are larger by 0.1 Å than the corresponding DD distances, which is in qualitative agreement with the experimental data, while the TTM3-F model predicts identical (within the statistical error) distances. By contrast, both sets of calculated OH and OD radial distribution functions show no shortening of the intramolecular OD bond with respect to OH, which are found to be 0.97 and 0.98 Å for the TTM2.1-F and TTM3-F models, respectively. Contrary to the experimental results, both water models do not also predict the lengthening of the H-bond in H2O (1.86 Å for both TTM2.1-F and TTM3-F) with respect to the D-bond in D2O (1.86 Å for TTM2.1-F and 1.85 Å for TTM3-F). Interestingly, shorter OD bonds and longer Hbonds have been predicted so far only by the quantum simulations within Car-Parrinello molecular dynamics of ref 154, which, however, have been shown to provide a very poor agreement with the measured ∆S(Q).63 At the same time, other quantum simulations with empirical force fields144,145 display similar trends as those predicted by the TTM2.1-F and TTM3-F models. 4. Dynamical Properties of Liquid Water

Figure 3. Comparison between the oxygen-oxygen radial distribution functions for liquid H2O (red) and D2O (black) computed at three different temperatures. Uncertainties are smaller than the thickness of the lines.

The dynamics of the H-bond network plays a central role in many fundamental phenomena occurring in Nature.35 In this regard, the behavior of light and heavy water may be very different. In order to characterize these differences at a molecular level, diffusive and relaxation processes in both liquids were extensively analyzed from ACMD simulations at T ) 298.15 K and T ) 278.15 K within the TTM2.1-F model.163 The self-diffusion coefficient was calculated using the wellknown expression169

D)

1 3

∫0∞ dt 〈v(t) · v(0)〉 31 ∫0∞ dt CVV(t)

(21)

where CVV is the velocity autocorrelation function. Insight into the relaxation processes at a molecular level was gained from the examination of the following autocorrelation functions

Cl(t) ) 〈Pl[e(t) · e(0)]〉

Figure 4. OH and HH radial distribution functions for liquid H2O (left panels) and D2O (right panels) computed at T ) 298.15 K with the TTM2.1-F (red) and TTM3-F (blue) models. Also shown are the experimental data from ref 38 (a, dashed black line) and from ref 64 (b, solid black line). Uncertainties in the computed results are smaller than the thickness of the lines.

corresponding to H-bonded (D-bonded) configurations at R ∼ 1.8 Å is shifted toward larger distances and its height is higher than the experimental estimates. Both the TTM2.1-F and TTM3-F models predict a symmetric shape for this peak in agreement with the data for H2O of ref 38. However, the most recent experimental results of ref 64 indicate that this prediction is only valid for D2O, while the measurements for light water resulted in an asymmetric distribution of

(22)

where Pl is the Legendre polynomial of order l ) 1, 2, 3 and e is a unit vector along the OH (OD) axis in the body-fixed reference frame of each water molecule. From integration of Cl(t), it is then possible to obtain the corresponding relaxation times, τl, which in the case of l ) 2 are directly accessible by experimental measurements. The quantum and classical results for the dynamical properties considered in this study are compared in Table 1 with the corresponding experimental data. First, in all cases, the quantum results are in better agreement with the experimental values, which stresses the fact that explicit inclusion of the effects arising from the quantization of the nuclear motion is necessary for a realistic description of the liquid properties when the underlying interactions are described by ab initio-based force fields, such as the TTM2.1-F model. However, it is important to note that the calculation of the diffusion coefficient using periodic boundary conditions is particularly sensitive to finite-



TABLE 1: Classical (C) and Quantum (Q) Results Compared with the Available Experimental Data for the Diffusion Coefficient (in Å2 ps-1) and for the Orientational Relaxation Times (in ps) for Liquid H2O and D2O at T ) 278.15 K and T ) 298.15 K (The Numbers in Parentheses Indicate the Statistical Error in Units of the Last Digit) H2O

D2O

T

exptl

C

Q

278.15

a

0.085(5) 10.05(5) 4.95(5) 3.40(5) 0.150(5) 6.45(5) 3.20(5) 2.15(5)

0.140(5) 6.25(5) 3.15(5) 2.10(5) 0.225(5) 4.70(5) 2.35(5) 1.55(5)

298.15

D τ1 τ2 τ3 D τ1 τ2 τ3

0.12

0.229b 2-7.5c 1.7-2.6d

exptl

0.19

H2O/D2O

C

Q

C

Q

0.080(5) 11.25(5) 5.45(5) 3.60(5) 0.140(5) 7.10(5) 3.45(5) 2.40(5)

0.110(5) 9.10(5) 4.50(5) 3.10(5) 0.185(5) 6.90(5) 3.45(5) 2.40(5)

1.1(1) 0.9(1) 0.9(1) 0.9(1) 1.1(1) 0.9(1) 0.9(1) 0.9(1)

1.3(1) 0.7(1) 0.7(1) 0.7(1) 1.2(1) 0.7(1) 0.7(1) 0.7(1)

a Extrapolated from ref 217 with statistical uncertainty equal to 5%. b From ref 218 with statistical uncertainty equal to 0.2%. c From refs 219-221. d From refs 192-194, 222.

size effects.170 In this regard, a value of D extrapolated to an infinite system can be derived from the following hydrodynamic relationship

D(L) ) D(∞) - ξ

kBT 6πηL

(23)

where η is the shear viscosity, L is the length of the simulation box, and ξ is a numerical coefficient that depends on the specific box geometry. In general, for a system consisting of 216 water molecules, D(∞) is about 10-20% larger than D(L), which means that, in the present case, the extrapolated quantum result will lie slightly above the experimental value. Second, explicit inclusion of quantum effects results in a faster water dynamics. For example, the self-diffusion coefficient of H2O increases by about 50% with respect to the classical result, which is consistent with the values reported in the literature for common empirical water force fields.15,134,145,146,171 By contrast, quantization of the D2O motion results only in a marginal increase of its self-diffusion coefficient, which indicates a more classical behavior of heavy water as expected from its relatively larger mass. Similar conclusions can also be reached from the analysis of the relaxation times. It is important to note that all previous quantum simulations with empirical force fields significantly overestimated the values of D and τl,145-147,171,172 which highlights again the importance of using ab initio-based force fields for a realistic and quantitative assessment of nuclear quantum effects in liquid water. Third, in line with the temperature dependence of the structural properties discussed in the previous section, the ratio between classical and quantum results for both H2O and D2O increases as the temperature decreases, which reinforces the notion that nuclear quantum effects become increasingly important at low temperatures. As will be shown in section 7, this may be particularly relevant for a realistic modeling of the ice properties. Importantly, the effects of the isotopic substitution on both the self-diffusion coefficient and relaxation times do not display any significant temperature dependence, which indirectly suggests that the contribution of tunneling, if present, should be small. 5. Infrared Spectrum and Hydrogen-Bond Dynamics in Liquid Water Infrared spectroscopy is an extremely powerful tool to characterize the structure and dynamics of liquid water. In Figure 5, the IR spectrum computed from ACMD simulations with the TTM3-F water model is compared with the experimental

Figure 5. Comparison between the calculated (red line) and the experimental (black line) infrared spectrum from ref 173 of liquid H2O at T ) 298.15 K.

measurements of ref 173. The calculated spectral features corresponding to the HOH bending (1550-1600 cm-1) and OH stretching (3000-3700 cm-1) vibrations appear to be slightly shifted by ∼100 cm-1 toward the red. However, their shapes are remarkably well reproduced. Somewhat less satisfactory is the agreement with the experimental band below 1000 cm-1 corresponding to the hindered motion of the water molecules. Clearly, the simulation underestimates the width of this broadband, which is directly connected to the differences between the calculated and experimental OH radial distribution functions in Figure 4. In addition, the shoulder at ∼200 cm-1 in the experimental IR spectrum, which has been shown to result from a dipole-induced dipole interaction,174 is also not accurately reproduced. We note here that the infrared spectrum for the TTM3-F water model was also computed using ring-polymer molecular dynamics, which however failed in reproducing the correct shape of the OH stretch band due to unphysical behavior of the ring-polymer MD method.138 A different development of the original TTM2.1-F model has also been proposed and applied with success to the investigation of the proton momentum distribution in water and ice.26 Since the OH frequencies of a water molecule are very sensitive to the surrounding environment, insight into the structure and dynamics of the liquid phase can be obtained from a detailed analysis of both position and shape of the corresponding IR band between 3000 and 3700 cm-1. From measurements of ultrafast IR spectroscopy, it was also suggested that tunneling can affect the reorientational dynamics of dilute

5710 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Figure 6. OH frequency distribution, P(ω01), computed for HOD in D2O. The three frequency regions shown in color correspond, respectively, to HOD molecules involved in strong (red), intermediate (green), and weak (blue) hydrogen bonds. See the main text for details.

HOD molecules in D2O.55 However, no frequency dependence of the anisotropy decay was found in similar IR pump-probe spectra of HOD in D2O with 45 fs pulses, which had enough bandwidth to span the entire OH absorption line.56 Although a great deal of information has been obtained about the dynamics of dilute HOD molecules in either H2O or D2O from different models based on classical mechanics,58,175-187 all of these simulations did not account for the intrinsic quantum motion of the water molecules. In order to fill this gap and to investigate in more detail the importance of quantum effects on the dynamics of the H-bond network, we have performed ACMD simulations at T ) 298.15 K for a system consisting of a single HOD molecule in a box of 127 D2O molecules at the experimentaldensityofD2O.Allsolute-solventandsolvent-solvent interactions (i.e., HOD-D2O and D2O-D2O interactions) were described by the TTM3-F water model, which, as illustrated in Figure 5 and discussed in more detail in ref 25, provides a realistic description of the IR band corresponding to the OH stretch as well as of other properties of liquid water. In total, 25 independent ACMD trajectories of 60 ps each were analyzed. Following ref 183, the instantaneous OH frequency of the HOD molecule was computed every 0.5 fs by solving numerically the Schro¨dinger equation for the corresponding one-dimensional OH oscillator. More specifically, for a given configuration of the HOD:D2O system, the OH potential energy curve, as defined by the TTM3-F model, was computed by stretching the OH bond while keeping the positions of all other atoms fixed. The second Legendre polynomial orientational correlation function, C2(t) in eq 22, was then computed in order to make a direct comparison with the results of pump-probe experiments. For this purpose, the distribution of the computed OH frequencies was divided into three subensembles corresponding to HOD molecules involved in H-bonds of different strengths (i.e., low, intermediate, and high frequency regions corresponding to strong, intermediate, and weak H-bonds). The OH frequency distribution, P(ω01), along with the three different frequency regions is shown in Figure 6. We note that the frequency distribution obtained from the quantum ACMD simulations with the TTM3-F model appears to be significantly different from analogous distributions computed from a combination of classical molecular dynamics simulations with empirical force fields and ab initio calculations (e.g., see ref 188). These differences will be analyzed in more detail in a forthcoming publication.

Paesani and Voth

Figure 7. Second Legendre polynomial orientational correlation function for the three different ranges of pump and probe OH frequencies defined in Figure 6. The dashed lines in black represent the corresponding single exponential fits. See the main text for details.

Three C2(t) (one for each of the frequency regions shown in Figure 6) were then computed for HOD molecules whose OH frequency lied in the same region at time 0 and time t, and the results are shown in Figure 7 using the same color scheme employed in Figure 6. All three C2(t) show an initial drop, which corresponds to the inertial motion of the HOD molecule and, consequently, depends on the initial OH frequency, as discussed in detail in ref 189. A weak frequency dependence is found up to t ∼ 1.0 ps, after which all three curves display similar longtime decays that can effectively be interpolated by a single exponential (black dashed lines in the figure) with a relaxation time of τ2 ) 2.0 ps-1. The present quantum results thus suggest that no tunneling is at play in the reorientation dynamics of HOD in D2O at ambient conditions, which is consistent with the molecular jump mechanism for hydrogen-bond breaking and forming.190 The frequency dependence displayed by C2(t) in the first picosecond after the excitation is similar to that observed in the pump-probe experiments of ref 54 and also found recently in analogous spectroscopic measurements for the HOD orientational dynamics in D2O.57 The behavior of C2(t) within the first picosecond suggests that further theoretical investigations on the relationship between H-bond configurations and orientational dynamics are needed. Importantly, the calculated value of τ2 is slightly smaller than the experimental estimates reported in the literature191-194 as well as the relaxation time computed with the TTM2.1-F model (see Table 1), which suggests that the water dynamics predicted by the TTM3-F model is probably too fast. In this regard, we note that the selfdiffusion coefficient computed from quantum simulations with the TTM3-F model25 also appears to be significantly larger than that previously obtained with the TTM2.1-F model.163 Given the recent debate over the average three-dimensional arrangement of the water molecules, we also investigated the occurrence of different topologies of H-bond configurations during the ACMD simulations from an analysis of their contributions to the OH stretch band. For this purpose, we followed the approach proposed in ref 188. First, the OH frequency distribution of Figure 6 was decomposed on the basis of the possible hydrogen-bond configurations of the HOD molecule, with the H-bond configurations being determined from the H-bond definition of ref 40. We then described each H-bond



TABLE 2: Analysis of the Different Topologies (Classes) of H-Bond Configurations Involving the HOD Molecule in the ACDM Simulations of the HOD:D2O System with the TTM3-F Water Modela NO

NH

ND

label

〈ω01〉

f

1 1 1 1 2 2 2 2

0 0 1 1 0 0 1 1

0 1 0 1 0 1 0 1

1N 2SD 2SH 3D 2N 3SD 3SH 4D

3431 3445 3294 3311 3383 3384 3246 3253

0.024 0.078 0.074 0.228 0.017 0.087 0.085 0.407

a See the main text for a description of NO, NH, and ND, and associated labels. 〈ω01〉 is the average frequency (in cm-1) of each class, and f is the fraction of HOD molecules in each class.

Figure 8. Distribution of the OH stretch frequencies computed from ACMD simulations of the HOD:D2O systems at T ) 298.15 K for the eight H-bond classes involving the HOD molecule.

configuration on the basis of the number of H-bonds involving the O atom (nO), the H atom (nH), and the D atoms (nD) of the HOD molecule. Eight classes of H-bond configurations were then defined according to the triplet of numbers NO, NH, and ND. The possible values are NO ) 1 corresponding to nO ) 0, 1, and NO ) 2 corresponding to nO ) 2, 3; NH ) 0 corresponding to nH ) 0, and NH ) 1 corresponding to nH ) 1, 2; ND ) 0 corresponding to nD ) 0, and ND ) 1 corresponding to nD ) 1, 2. As described in ref 188, these classes were labeled according to the total number of H-bonds, the number of H-bond donors, and, the H or D label in the case of single donors. Table 2 lists all eight classes along with the corresponding labels and the fraction of configurations in each class, while Figure 8 shows the frequency distribution for each class, together with the overall distribution of Figure 6. The present results are in qualitative agreement with those from the classical simulations of ref 188 and, in particular, share with them the same order of the OH frequencies as a function of the number of H-bonds. However, significant differences exist for the shape of each frequency distribution, which is most likely due to the inclusion of nuclear quantum effects in the present calculations as well as to the different schemes employed to compute the instantaneous OH frequencies. Nevertheless, on the basis of the data of Table 2 and Figure 8, similar conclusions to those of ref 188 can be drawn. More specifically, the present ACMD simulations with the TTM3-F model indicate that the fraction of water molecules with only one donating H-bond is smaller than that

of double donors, supporting the usual picture of liquid water in which each molecule is, on average, arranged in a tetrahedral geometry. The present quantum analysis also indicates that some refinements of the TTM water models are likely necessary in order to obtain a more accurate description of the properties of liquid water. In particular, it appears that both the TTM2.1-F and TTM3-F models predict H-bonds that are too long, which may affect the dynamics of the hydrogen-bond network as well as some thermodynamic properties (e.g., the enthalpy of vaporization), which have also been shown to not be well reproduced.25,163,195 6. Solvent Isotope Effects on the Stability and Hydration of Biological Systems As mentioned in the Introduction, H2O and D2O may have very different effects on the properties of biological systems. Since it has been shown that protein motion is intimately linked to the motion of the surrounding aqueous environment,196-198 a molecular-level characterization of solvent isotope effects may help in developing a more comprehensive understanding of the relationship between protein structure and function under physiologically relevant conditions. Furthermore, understanding solvent isotope effects on the properties of biomolecules is also of importance from a practical point of view. For example, liquid D2O is the solvent of choice in many experiments aimed at determining mechanistic aspects of protein-folding kinetics. Pulsed hydrogen/deuterium exchange (HDX) methods coupled to mass spectrometry indeed exploit the possibility of the hydrogen atoms of OH, SH, and NH groups of a protein to exchange with the deuterium atoms of the solvent.199 Since H/D exchanges are only possible if the H atoms are available to the solvent and not involved in stable hydrogen bonds, monitoring the rate of these exchanges in time provides a route to identify the most likely folding pathways. Moreover, since the substitution of H with D also affects the thermodynamic and kinetic constants, liquid D2O is commonly used to elucidate the mechanism of enzymatic reactions.200 From a theoretical point of view, very little is known about the effects of D2O on the properties of biomolecules, and previous studies were limited to a modeling of the isotope effects based on macroscopic thermodynamics.201 The reason for this may have its origins in the lack of an accurate and efficient theoretical and computational approach that could allow for a quantitative investigation at a molecular level of nuclear quantum effects in realistic biological systems. In order to address this issue, we have recently developed and implemented a completely new module in the AMBER suite of codes for biomolecular simulations,202 which allows the user to perform large-scale quantum simulations of complex biological systems within the path-integral molecular dynamics and centroid molecular dynamics formalisms. As a first step toward a quantitative assessment of the solvent isotope effect on the properties of biomolecules, the stability of five amino acids (AAs), whose physicochemical properties are representative of those of the different naturally occurring amino acids, and their corresponding hydration structures were investigated in both liquid H2O and D2O. More specifically, we analyzed the behavior of two hydrophobic amino acidssalanine (aliphatic) and phenylalanine (aromatic)sand three polar amino acidssserine (neutral), glutamate (negatively charged), and arginine (positively charged). All simulations were carried out for a system consisting of a single amino acid and 256 water molecules at T ) 298.15 K using the NMPIMD implementation


Paesani and Voth

in version 10 of AMBER.202 The interactions between the water molecules were described by the quantum q-SPC/Fw force field,15 while the GAFF force field was used for all amino acids.203 In the latter case, the partial charges were obtained following the standard RESP procedure,204 and each amino acid was terminated by the standard NME and ACE protecting groups. In all cases, after 2 ns of equilibration, the simulation was run for another 10 ns, the data from which was used for the actual calculation of both thermodynamic and structural properties. The relative stability of the five amino acids in heavy and light water was analyzed by computing the free energy change (∆∆F) associated with their transfer from H2O to D2O. This free energy change, which is proportional to the ratio of the corresponding quantum partition functions, can be expressed as

(

Q(AA+D2O) Q(AA+H2O) Q(D2O) -kT ln (H O) (24) Q 2

AA+D2O ∆∆F ) ∆FAA+H - ∆FHD22OO ) -kT ln 2O

(

)

)

where the first term corresponds to the free energy difference between the two combined systems AA:D2O and AA:H2O, while the second term corresponds to the free energy difference between pure D2O and H2O. Each term was computed using the so-called thermodynamic integration with respect to mass,205 which is formulated within the path-integral formalism as

∆F )

∫0

1

〈

Figure 9. Free energy of transfer (in kcal mol-1) from H2O to D2O computed according to eq 24.

〉

dVeff(λ) dλ dλ

(25)

where Veff(λ) is given by

Veff(λ) ) -kT ln Q(λ)

(26)

and λ is a parameter that interpolates between the mass of the hydrogen (λ ) 0) and deuterium (λ ) 1) atoms

m(λ) ) (1 - λ)m(H) + λm(D)

(27)

Eleven values of λ with an interval of ∆λ ) 0.1 were employed to compute the integral in eq 25, and since we are specifically interested in the effects on the amino-acid properties from the isotopic substitution on the water molecules, only the H atoms of H2O were substituted with D atoms. In order to determine ∆∆F, the free energy difference between liquid H2O and D2O (second term in eq 24) was first calculated at T ) 298.15 K, resulting in a value of -4.15 ( 0.05 kcal mol-1. This is in remarkably good agreement with the corresponding experimental value of -4.13 kcal mol-1,206 which reinforces the notion that the q-SPC/Fw water model provides a reasonable description of the water properties at ambient conditions when explicit quantum effects are included in the simulation.15 Subsequently, the free energy of transfer from H2O to D2O was computed for all five amino acids according to eq 24. The results shown in Figure 9 indicate that ∆∆F is negative in all cases, displaying similar values that range between -0.02 ( 0.1 kcal mol-1 (alanine) and -0.35 ( 0.1 kcal mol-1 (glutamate). Qualitatively, the general trend appears to be related to the size

Figure 10. Structure of glutamate and alanine with atom labeling.

of the amino acids with the larger ones displaying a more negative value of ∆∆F. However, the statistical uncertainties as well as the small absolute values prevent us for a more quantitative assessment. Interestingly, the values of ∆∆F obtained for small amino acids (alanine and serine) are similar to those calculated for small hydrophobes using information theory.201 In order to investigate possible differences in the hydration structure, the radial distribution functions describing correlations between the oxygen atoms of the water molecules and the heavy atoms of glutamate and alanine were analyzed (see Figure 10 for the specific labeling of each heavy atom). Several RDFs computed for glutamate in either H2O or D2O are shown in Figure 11. Although both sets of results are qualitatively similar, there are still some differences, which are more evident in the RDFs involving the C3 and C4 carbon atoms of the side chain (see Figure 10). As a comparison, Figure 12 shows the RDFs computed for alanine. In this case, the differences between H2O and D2O are negligible and are barely visible in the RDFs involving C1 (see Figure 10). A more extensive analysis is clearly necessary to assess in a quantitative way the effects of H/D substitution on the behavior of realistic biological systems. In this regard, the sensitivity of the free energy calculations on the choice of the force field employed for the amino acids as well as a broader analysis of the solvent isotope effects on the hydration of organic and inorganic molecules will be the subject of a future study. It will also be important to investigate the temperature dependence of the free energy of transfer from H2O to D2O, which will possibly provide new insight into the role played by nuclear quantum effects on the “biological properties” of water. However, in light of the results discussed in the previous section regarding the nonlinear dependence of quantum effects with respect to the inverse temperature, we can anticipate that the q-SPC/Fw model, which was parametrized at T ) 298.15 K, may not be completely adequate for such studies and a more sophisticated (and



Figure 11. Radial distribution functions describing the correlation between the oxygen atoms of the water molecules and the heavy atoms of glutamate labeled as in Figure 10. H2O (red) and D2O (black). Uncertainties are smaller than the thickness of the lines.

Figure 12. Radial distribution functions describing the correlation between the oxygen atoms of the water molecules and the heavy atoms of alanine labeled as in Figure 10. H2O (red) and D2O (blue). Uncertainties are smaller than the thickness of the lines.

transferable) water force field (e.g., the TTM water models) might be necessary. 7. Structural and Dynamical Properties of the Ice Surface The ice surface plays a central role in many physical chemical processes, some of which are especially relevant to the variations of the atmosphere composition and, consequently, to the Earth’s climate.207 For this reason, several theoretical studies that

combined classical molecular dynamics with empirical water force fields, as well as ab initio molecular dynamics simulations, have been performed to characterize at a molecular level the properties of the ice surface for T < 273.15 K.208-212 In particular, extensive simulations of the (0001) surface with the TIP4P water model showed that a disordered phase in the topmost layers appears for T > 230 K.208 However, it has recently been demonstrated that the actual melting point for the TIP4P water model corresponds to T0 ) 232 K. After a proper scaling of


Figure 13. Quantum and classical oxygen-oxygen radial distribution functions computed for the solid-liquid coexistence of water as a function of the temperature. In red are the RDFs corresponding to molecules initially located in the ice phase, while in blue are those for molecules initially located in the water phase. Uncertainties are smaller than the thickness of the lines.

the temperatures with respect to the actual melting point of the TIP4P model, it is possible to see that the simulations of ref 208 do not predict the appearance of a liquid-like layer at the ice surface. More recently, using a different definition of the liquid-like behavior for the water molecules, classical simulations have shown the appearance of a thin QQL below the actual melting point computed for several empirical water models.213 In all of these simulations, nuclear quantum effects were neglected, which, however, should become more and more significant as the temperature decreases. Therefore, inclusion of such effects in simulations of ice below the melting point is particularly important for a realistic modeling of its properties. Recently, we have carried out an extensive analysis of the temperature dependence of the thermodynamic and dynamical properties of the ice surface with explicit inclusion of the effects arising from the quantization of the nuclear motion.214 The quantum simulations were performed with NMPIMD and CMD methods, while the interactions between water molecules were described by the ab initio-based, polarizable TTM2.1-F model. For comparison, we also performed standard classical molecular dynamics simulations. The actual melting point of ice Ih modeled with the TTM2.1-F water force field was calculated from simulations of the solid-liquid coexistence for a system consisting of 192 water molecules that were evenly distributed in the two phases at the beginning of both NMPIMD and classical simulations. The oxygen-oxygen radial distribution function was monitored as the temperature was raised from T ) 200 K in increments of ∆T ) 10 K (see Figure 13). After an initial rough estimate of T0, smaller increments of 2.5 K were used to determine the melting temperature with better precision. The quantum NMPIMD simulations provided a melting point of T0 ) 227.5 ( 1.5 K, while T0 ) 242.5 ( 1.5 K was estimated from the corresponding classical simulations. The difference from the experimental value (T0 ) 273.15 K) arises from a slight underestimation of the enthalpy of fusion. For quantum ice, NMPIMD simulations indeed provided a value of ∆Hf ) 1.0 ( 0.1 kcal mol-1, while ∆Hf ) 1.1 ( 0.1 kcal mol-1 was obtained from the corresponding classical calculations, which are smaller than the experimental value of ∆Hf ) 1.44 kcal mol-1.215

Paesani and Voth

Figure 14. Top (a) and side (b) views of the model of hexagonal ice employed in all simulations of the ice surface. Large red spheres represent oxygen atoms, while small white spheres correspond to the hydrogen atoms.

The investigation of the properties of the ice surface below the melting point were then carried out as a function of the relative temperature difference from T - T0 ) -77.5 K to T T0 ) -2.5 K for a slab geometry representing the (0001) surface, which consisted of 192 water molecules placed in 16 monolayers (Figure 14). In order to simulate a semi-infinite system, the two bottom-most layers were kept fixed throughout the simulations. The basal (0001) surface was chosen because it is considered to be the most prevalent ice surface under polar stratospheric conditions and, consequently, is thought to play a key role in the catalysis of chlorine activation reactions.91 Two order parameters were calculated to quantify the degree of surface disorder as well as the influence of nuclear quantum effects on the properties of the ice surface. First, a translational order parameter was computed as 3

ST )

12

12

∑ ∑ ∑ cos(km · Rij)/NL2

1 3 m)1

(28)

i)1 j)1

where k1 ) (2π/a)(1,1/√3,0), k2 and k3 are related to k1 by one or two subsequent rotations of 120°,208 and NL ) 12 is the number of water molecules initially located in each layer. In eq 28, Rij is a vector pointing from the oxygen atom of molecule i to the oxygen atom of molecule j belonging to the same layer, and a is the lattice constant that defines the (0001) plane. For a completely ordered layer, ST ) 1, while, for a disordered system, ST ) 1/NL ) 1/12 ) 0.084. Second, the rotational order parameter NL

SR )

∑ 2.5(xi4 + yi4 + zi4 - 0.6)/NL

(29)

i)1

was computed following the same procedure described in ref 212. The initial dipole vector of molecule i was normalized at the beginning of the simulation, and a rotation matrix was obtained by transforming the normalized vector to a unit vector parallel to the z-axis. During the simulation, the dipole vector of molecule i was then normalized and multiplied by the same rotation matrix to obtain xi, yi, and zi. For an ordered lattice, SR ) 1, while SR f 0 for a completely disordered system.



TABLE 3: Translational Order Parameter (ST) Computed from Quantum (First Entry) and Classical (Second Entry) Simulations of the Ice Surface As a Function of the Relative Temperature Difference T - T0

TABLE 4: Rotational Order Parameter (SR) Computed from Quantum (First Entry) and Classical (Second Entry) Simulations of the Ice Surface As a Function of the Relative Temperature Difference T - T0

T - T0 Nlayer 3 4 5 6 7 8 9 10 11 12 13 14 15 16

T - T0

-77.5 K

-57.5 K

-37.5 K

-27.5 K

-17.5 K

-7.5 K

-2.5 K

0.93 0.94 0.93 0.94 0.92 0.93 0.92 0.93 0.92 0.93 0.92 0.93 0.92 0.93 0.92 0.93 0.91 0.92 0.91 0.92 0.91 0.92 0.91 0.92 0.88 0.90 0.88 0.89

0.93 0.93 0.93 0.93 0.91 0.92 0.91 0.92 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.90 0.91 0.90 0.91 0.90 0.91 0.88 0.90 0.85 0.87 0.83 0.84

0.91 0.92 0.91 0.92 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.89 0.89 0.89 0.89 0.89 0.88 0.89 0.87 0.88 0.85 0.84 0.80 0.82

0.91 0.92 0.91 0.92 0.89 0.90 0.89 0.90 0.89 0.90 0.89 0.90 0.89 0.90 0.89 0.90 0.87 0.89 0.87 0.89 0.86 0.89 0.85 0.89 0.24 0.37 0.22 0.35

0.90 0.91 0.90 0.91 0.89 0.90 0.89 0.89 0.88 0.89 0.88 0.89 0.88 0.89 0.88 0.89 0.87 0.88 0.87 0.88 0.85 0.88 0.84 0.88 0.10 0.38 0.09 0.34

0.90 0.91 0.90 0.91 0.88 0.89 0.88 0.89 0.88 0.88 0.87 0.88 0.87 0.88 0.87 0.88 0.87 0.88 0.87 0.88 0.83 0.81 0.80 0.65 0.10 0.12 0.09 0.10

0.90 0.91 0.90 0.91 0.88 0.89 0.88 0.89 0.87 0.88 0.87 0.88 0.87 0.88 0.87 0.88 0.68 0.87 0.64 0.87 0.10 0.81 0.10 0.64 0.09 0.10 0.09 0.09

The quantum and classical results for ST and SR are reported in Tables 3 and 4, respectively.214 Analysis of these two order parameters as a function of the temperature and layer location along the z-direction indicates that an ordered structure is essentially conserved at all temperatures for the most inner layers (Nlayer e 12). However, it is interesting to note that the quantum simulations already predict some degree of disorder for layers 11 and 12 at the highest temperatures. Qualitative differences between the quantum and classical results are instead evident for the four topmost layers as the temperature increases above T - T0 > -30 K. In particular, a translational and rotational disorder appears in all four topmost layers of quantum ice at T - T0 ) -2.5 K. By contrast, only the two topmost layers show a similar disorder in the classical case. More importantly, disorder in the two topmost layers of quantum ice persists down to a temperature of T - T0 ) -27.5 K, which is about 20 K lower than the limit obtained from the corresponding classical simulations. The difference between the quantum and classical models of the ice surface at T - T0 ) -2.5 K can be clearly recognized from the inspection of the density profiles shown in Figure 15.214 These curves were obtained by computing the distance of each oxygen atom from the bottom-most layer along a direction perpendicular to the (0001) plane (the z-direction). The quantum results predict the existence of a liquid-like structure at the nearsurface region that comprises four water layers and extends for more than 10 Å. This value is in fairly good agreement with estimates of the QLL thickness estimated in optical ellipsometry experiments216 and photoelectron spectroscopy.217 In line with the results for SR and ST, the classical profile instead displays a structural disorder only in the two topmost layers.

Nlayer 3 4 5 6 7 8 9 10 11 12 13 14 15 16

-77.5 K

-57.5 K

-37.5 K

-27.5 K

-17.5 K

-7.5 K

-2.5 K

0.80 0.77 0.76 0.77 0.72 0.77 0.72 0.77 0.72 0.77 0.72 0.77 0.72 0.77 0.72 0.77 0.72 0.77 0.72 0.75 0.71 0.74 0.71 0.74 0.65 0.66 0.53 0.56

0.77 0.75 0.72 0.73 0.70 0.72 0.70 0.72 0.69 0.72 0.69 0.72 0.69 0.72 0.69 0.72 0.69 0.72 0.68 0.72 0.67 0.70 0.66 0.66 0.61 0.61 0.46 0.51

0.77 0.77 0.71 0.72 0.68 0.71 0.67 0.70 0.66 0.70 0.65 0.70 0.65 0.70 0.65 0.70 0.65 0.69 0.65 0.68 0.64 0.67 0.61 0.62 0.55 0.54 0.40 0.42

0.71 0.75 0.66 0.72 0.64 0.69 0.62 0.68 0.61 0.66 0.60 0.66 0.60 0.66 0.60 0.65 0.60 0.67 0.60 0.65 0.57 0.65 0.50 0.60 0.06 0.42 0.04 0.31

0.69 0.73 0.66 0.71 0.64 0.70 0.62 0.69 0.61 0.68 0.60 0.68 0.60 0.67 0.60 0.66 0.60 0.66 0.60 0.63 0.59 0.60 0.52 0.58 0.06 0.40 0.04 0.25

0.69 0.73 0.65 0.70 0.63 0.68 0.62 0.66 0.61 0.65 0.61 0.65 0.61 0.64 0.58 0.62 0.58 0.61 0.57 0.58 0.52 0.55 0.42 0.44 0.04 0.04 0.03 0.03

0.66 0.69 0.64 0.67 0.62 0.66 0.60 0.66 0.60 0.65 0.58 0.64 0.57 0.63 0.55 0.62 0.45 0.60 0.34 0.58 0.09 0.55 0.07 0.35 0.05 0.04 0.04 0.03

This quasi-liquid structure of the topmost water layers is also accompanied by a larger lateral mobility of H2O molecules, which can be quantified by the diffusion coefficient (D) for motions parallel to the surface. Quantum simulations,214 which were carried out using the ACMD method, provided a value of D| ) 2.1 Å2 ps-1 which is significantly larger than the corresponding classical value of D| ) 1.2 Å2 ps-1. It is important to note that the classical results214 are similar to those obtained in ref 208 after scaling the temperature with respect to the actual melting point of the TIP4P model. This suggests that the qualitatively different behavior seen for the

Figure 15. Quantum (red) and classical (blue) density profiles along the z-direction perpendicular to the (0001) ice surface computed at T - T0 ) -2.5 K. Uncertainties are smaller than the thickness of the lines.

5716 J. Phys. Chem. B, Vol. 113, No. 17, 2009 structural and dynamical properties of the ice surface computed at a quantum and a classical level should be mostly independent of the water model employed in the simulations and effectively due to the explicit inclusion of the effects of quantization of the nuclear motion. Such differences may have important implications in the characterization of the physicochemical processes occurring at the ice surface. In particular, the prediction from quantum simulations of a thicker quasi-liquid layer may be a key factor to understand the role played by the ice surface of polar stratospheric clouds in the catalysis of chlorine activation reactions, which are suggested to lead to the observed ozone depletion in the stratosphere. 8. Conclusions In this review, we have discussed recent progress toward the development of an efficient and reliable approach for investigating the properties of water at a molecular level under different conditions and in different environments. In particular, we have focused our attention on theoretical and computational methodologies that combine a proper treatment of the nuclear dynamics at a quantum-mechanical level using path-integral molecular dynamics and centroid molecular dynamics, with an accurate description of the underlying Born-Oppenheimer potential energy surface. We have analyzed the structural and dynamical properties of liquid H2O and D2O as a function of the temperature from T ) 353.15 K to T ) 273.15 K. A comparison of the quantum and classical results for the intensity of X-ray scattering shows that the quantum simulations provide very good agreement with the corresponding experimental data over the entire range of temperatures, indicating that nuclear quantum effects become increasingly important for an accurate description of the liquid structure near the melting temperature. Inspection of the differences between the structure factors for liquid H2O and D2O computed at three different temperatures reveal a relatively larger probability for H2O molecules to be located between the first and second solvation shell, which was confirmed by a direct comparison of the corresponding oxygen-oxygen radial distribution functions. Analysis of the water dynamical properties also indicates that explicit inclusion of nuclear quantum effects is important to reproduce the experimental values. The absence of any significant temperature dependence in the ratio between the dynamical properties computed for D2O and H2O suggested a negligible role for hydrogen tunneling. This behavior was confirmed by a direct simulation of the dynamics of dilute HOD in D2O carried out for the first time at a quantum-mechanical level using a centroid molecular dynamics. It has also been shown that the quantum simulations with ab initio-based water models predict an ensemble of different local structures for the liquid with a prevalence of configurations in which each molecule is, on average, 4-fold coordinated. We have also investigated the relative stability of different amino acids and their corresponding hydration structure in light and heavy water at ambient conditions using path-integral molecular dynamics in combination with an empirical quantum model for water. The free energies of transfer from H2O to D2O calculated for each amino acid examined in this study lie between -0.02 kcal mol-1 for alanine and -0.35 kcal mol-1 for glutamate, and appear to be sensitive to the size of the amino acids. Inspection of the hydration structure around polar (glutamate) and hydrophobic (alanine) amino acids shows some appreciable differences between H2O and D2O, which may be relevant for a more precise interpretation of experimental investigations carried out in heavy water.

Paesani and Voth Finally, the calculated quantum structural and dynamical properties of the (0001) surface of ice Ih below the melting point have been reviewed.214 This analysis has shown that an explicit inclusion of nuclear quantum effects has a significant impact on the description of the ice surface, resulting in a picture quite different from the one obtained from purely classical simulations. In particular, the quantum results obtained at 2.5 K below the simulated melting point predict the existence of a much thicker quasi-liquid phase, which comprises four ice layers and extends for approximately 10 Å into the bulk, in good agreement with estimates from optical ellipsometry and photoelectron experiments. Importantly, translational and rotational disorder in the two topmost layers of quantum ice persists at temperatures that are about 10-20 K lower than their classical counterparts, which should have important implications in the modeling of the heterogeneous processes occurring on the ice surface. Although it has been shown that explicit inclusion of nuclear quantum effects in simulations with ab initio-based TTM water models generally improves the agreement with the available experimental data, there are still some discrepancies. In particular, the length of the hydrogen bonds appears to be slightly overestimated, which may have direct repercussions on the thermodynamic and dynamical properties. For example, the enthalpy of vaporization of liquid water is overestimated, while the enthalpy of fusion is underestimated, the latter resulting in a predicted value for the melting point of ice that is too low. In addition, the differences in the hydrogen-bond geometry may also be responsible for the discrepancies found in the infrared spectrum below 1000 cm-1, and in a somewhat too fast orientational dynamics of the HOD molecule in D2O. Some refinements of the TTM force field thus appear to be necessary in order to develop a universal water model capable of quantitatively describing water properties from clusters to bulk, and under different conditions of temperature and pressure. In the future, we plan to study the solvent isotope effects on the hydration structure of more realistic systems (e.g., small peptides) as well as to investigate more specifically the hydrophobic effect in either liquid H2O or D2O as a function of temperature. Moreover, the quantum approach described here to analyze the dynamics of the hydrogen-bond network in bulk will be extended to examine the behavior of water molecules at various interfaces with a particular emphasis on those of atmospheric and biological relevance. In this context, the development of general ab initio-based force fields, such as the TTM water model, will certainly be critical for achieving a quantitative, molecular-level description and understanding of these systems through quantum simulations. Acknowledgment. This research was supported by the Office of Naval Research through grant N00014-05-1-0457. The authors thank Prof. Dave Case and Dr. Wei Zhang for their help with the development and implementation of nuclear quantum methods in AMBER, Dr. Sotiris Xantheas for many discussions about the TTM water models, and Drs. Satoru Iuchi, VinodKrishna,andKimWongforvaluablescientificconversations. References and Notes (1) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (2) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (3) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269.

Centennial Feature Article (4) Reimers, J. R.; Watts, R. O.; Klein, M. L. Chem. Phys. 1982, 64, 95. (5) 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. (6) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (7) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123. (8) Abascal, J. L. F.; Sanz, E.; Fernandez, R. G.; Vega, C. J. Chem. Phys. 2005, 122. (9) Kumar, R.; Skinner, J. L. J. Phys. Chem. B 2008, 112, 8311. (10) Dang, L. X.; Pettitt, B. M. J. Phys. Chem. 1987, 91, 3349. (11) Ferguson, D. M. J. Comput. Chem. 1995, 16, 501. (12) Teleman, O.; Jonsson, B.; Engstrom, S. Mol. Phys. 1987, 60, 193. (13) Amira, S.; Spangberg, D.; Hermansson, K. Chem. Phys. 2004, 303, 327. (14) Wu, Y. J.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (15) Paesani, F.; Zhang, W.; Case, D. A.; Cheatham, T. E., III; Voth, G. A. J. Chem. Phys. 2006, 125, 184507. (16) Lobaugh, J.; Voth, G. A. J. Chem. Phys. 1997, 106, 2400. (17) Dang, L. X.; Chang, T. M. J. Chem. Phys. 1997, 106, 8149. (18) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (19) Saint-Martin, H.; Hernandez-Cobos, J.; Bernal-Uruchurtu, M. I.; Ortega-Blake, I.; Berendsen, H. J. C. J. Chem. Phys. 2000, 113, 10899. (20) Lamoureux, G.; MacKerell, A. D. J.; Roux, B. J. Chem. Phys. 2003, 119, 5185. (21) Burnham, C. J.; Li, J. C.; Xantheas, S. S.; Leslie, M. J. Chem. Phys. 1999, 110, 4566. (22) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1500. (23) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115. (24) Fanourgakis, G. S.; Xantheas, S. S. J. Phys. Chem. A 2006, 110, 4100. (25) Fanourgakis, G. S.; Xantheas, S. S. J. Chem. Phys. 2008, 128, 074506. (26) Burnham, C. J.; Anick, D. J.; Mankoo, P. K.; Reiter, G. F. J. Chem. Phys. 2008, 128, 154519. (27) Jeon, J.; Lefohn, A. E.; Voth, G. A. J. Chem. Phys. 2003, 118, 7504. (28) Iuchi, S.; Morita, A.; Kato, S. J. Phys. Chem. B 2002, 106, 3466. (29) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933. (30) Donchev, A. G.; Galkin, N. G.; Illarionov, A. A.; Khoruzhii, O. V.; Olevanov, M. A.; Ozrin, V. D.; Subbotin, M. V.; Tarasov, V. I. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 8613. (31) Bukowski, R.; Szalewicz, K.; Groenenboom, G. C.; van der Avoird, A. Science 2007, 315, 1249. (32) Goldman, N.; Leforestier, C.; Saykally, R. J. Philos. Trans. R. Soc., Ser. A 2005, 363, 493. (33) Mankoo, P. K.; Keyes, T. J. Chem. Phys. 2008, 129, 034504. (34) Dang, L. X. J. Phys. Chem. B 1998, 102, 620. (35) Marećhal, Y. The Hydrogen Bond and the Water Molecule: The Physics and Chemistry of Water, Aqueous and Bio-Media; Elsevier: Amsterdam, The Netherlands, 2006. (36) Kennedy, D.; Norman, C. Science 2005, 309, 75. (37) Sorenson, J. M.; Hura, G.; Glaeser, R. M.; Head-Gordon, T. J. Chem. Phys. 2000, 113, 9149. (38) Soper, A. K. Chem. Phys. 2000, 258, 121. (39) Modig, K.; Pfrommer, B. G.; Halle, B. Phys. ReV. Lett. 2003, 90, 075502. (40) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Na¨slund, L. Å.; Hirsh, T. K.; Ojamaë, L.; Glatzel, P.; Pettersson, L. G. M.; Nilsson, A. Science 2004, 304, 995. (41) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Science 2004, 306, 851. (42) Head-Gordon, T.; Johnson, M. E. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 7973. (43) Mantz, Y. A.; Chen, B.; Martyna, G. J. Chem. Phys. Lett. 2005, 405, 294. (44) Mantz, Y. A.; Chen, B.; Martyna, G. J. J. Phys. Chem. B 2006, 110, 3540. (45) Soper, A. K. J. Phys.: Condens. Matter 2005, 17, S3273. (46) Leetmaa, M.; Ljungberg, M.; Ogasawara, H.; Odelius, M.; Na¨slund, L.-A.; Nilsson, A.; Pettersson, L. G. M. J. Chem. Phys. 2006, 125, 244510. (47) Prendergast, D.; Galli, G. Phys. ReV. Lett. 2006, 96, 215502. (48) Leetmaa, M.; Wikfeldt, K. T.; Ljungberg, M. P.; Odelius, M.; Swenson, J.; Nilsson, A.; Pettersson, L. G. M. J. Chem. Phys. 2008, 129, 084502. (49) Tokushima, T.; Harada, Y.; Takahashi, O.; Senba, Y.; Ohashi, H.; Pettersson, L. G. M.; Nilsson, A.; Shin, S. Chem. Phys. Lett. 2008, 460, 387. (50) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Science 2004, 306, 851.

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5717 (51) Wang, Z. H.; Pakoulev, A.; Pang, Y.; Dlott, D. D. Chem. Phys. Lett. 2003, 378, 281. (52) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Cohen, R. C.; Geissler, P. L.; Saykally, R. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 14171. (53) Walrafen, G. E.; Fisher, M. R.; Hokmabadi, M. S.; Yang, W. H. J. Chem. Phys. 1986, 85, 6970. (54) Nienhuys, H. K.; van Santen, R. A.; Bakker, H. J. J. Chem. Phys. 2000, 112, 8487. (55) Rezus, Y. L. A.; Bakker, H. J. J. Chem. Phys. 2005, 123, 114502. (56) Loparo, J. J.; Fecko, C. J.; Eaves, J. D.; Roberts, S. T.; Tokmakoff, A. Phys. ReV. B 2004, 70. (57) Bakker, H. J.; Rezus, Y. L. A.; Timmer, R. L. A. J. Phys. Chem. A 2008, 112, 11523. (58) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 264. (59) Laage, D.; Hynes, J. T. Chem. Phys. Lett. 2006, 433, 80. (60) Laage, D.; Hynes, J. T. J. Phys. Chem. B 2008, 112, 14230. (61) Tauber, M. J.; Mathies, R. A. J. Am. Chem. Soc. 2003, 125, 1394. (62) Hart, R. T.; Benmore, C. J.; Neuefeind, J. C.; Kohara, S.; Tomberli, B.; Egelstaff, P. A. Phys. ReV. Lett. 2005, 94, 047801. (63) Hart, R. T.; Mei, Q.; Benmore, C. J.; Neuefeind, J. C.; Turner, J. F. C.; Dolgos, M.; Tomberli, B.; Egelstaff, P. A. J. Chem. Phys. 2006, 124, 134505. (64) Soper, A. K.; Benmore, C. J. Phys. ReV. Lett. 2008, 101, 065502. (65) Kushner, D. J.; Baker, A.; Dunstall, T. G. Can. J. Physiol. Pharmacol. 1999, 77, 79. (66) Kresheck, G. C.; Schneider, H.; Scheraga, H. A. J. Phys. Chem. 1965, 69, 3132. (67) Marcus, Y.; Ben-Naim, A. J. Chem. Phys. 1985, 83, 4744. (68) Maybury, R. H.; Katz, J. J. Nature 1956, 177, 629. (69) Hermans, J. J.; Sheraga, H. A. Biochim. Biophys. Acta 1959, 36, 534. (70) Dong, A.; Kendrick, B.; Kreigard, L.; Matsuura, J.; Manning, M. C.; Capenter, J. F. Arch. Biochem. Biophys. 1998, 391. (71) Parker, M. J.; Clarke, A. R. Biochemistry 1997, 36, 5786. (72) Kern, D.; Zaccai, G.; Giege, R. Biophys. Chem. 1980, 59, 203. (73) Makhatadze, G. I.; Clore, G. M.; Gronenborn, A. M. Nat. Struct. Biol. 1995, 2, 852. (74) Bonnete, F.; Madern, D.; Zaccai, G. J. Mol. Biol. 1994, 244, 436. (75) Baghurst, P. A.; Nichol, L. W.; Sawyer, W. H. J. Biol. Chem. 1972, 247, 3199. (76) Harmony, J. A.; Himes, R. H.; Schowen, R. L. Biochemistry 1975, 14, 5379. (77) Woodfin, B. M.; Henderson, R. F.; Henderson, T. R. J. Biol. Chem. 1970, 245, 3733. (78) Omori, H.; Kuroda, M.; Naora, H.; Takeda, H.; Nio, Y.; Otani, H.; Tamura, K. J. Cell Biol. 1997, 74, 273. (79) Panda, D.; Chakrabarti, G.; Hudson, J.; Pigg, K.; Miller, H. P.; Wilson, L.; Himes, R. H. Biochemistry 2000, 39, 5075. (80) Chakrabarti, G.; Kim, B.; Gupta, M. L.; Barton, J. S., Jr.; Himes, R. H. Biochemistry 1999, 38, 3067. (81) Cioni, P.; Strambini, G. B. Biophys. J. 2002, 82, 3246. (82) Ikeda, M.; Suzuki, S.; Kishio, M.; Hirono, M.; Sugiyama, T.; Matsuura, J.; Suzuki, T.; Sota, T.; Allen, C. N.; Konishi, S.; Yoshioka, T. Biophys. J. 2004, 86, 565. (83) Altermatt, H. J.; Gebbers, J. O.; Laissue, J. A. Cancer 1988, 51, 303. (84) Takeda, H.; Nio, Y.; Omori, H.; Uegaki, K.; Hirahara, N.; Sasaki, S.; Tamura, K.; Ohtani, H. Anti-Cancer Drugs 1998, 9, 715. (85) Bhak, J. Y.; Lee, J.-H.; Chung, H. S.; Lee, H. Y.; Chung, B. C.; Park, M. S.; Min, S. K.; Kim, M. O. J. Ind. Eng. Chem. 2007, 13, 501. (86) Dash, J. G.; Fu, H.; Wettlaufer, J. S. Rep. Prog. Phys. 1995, 58, 115. (87) Wettlaufer, J. S. Philos. Trans. R. Soc. London, Ser. A 1999, 357, 3403. (88) Li, Y.; Somorjai, G. A. J. Phys. Chem. C 2007, 111, 9631. (89) Rosenberg, R. Phys. Today 2005, 58, 50. (90) Makkonen, L. J. Phys. Chem. B 1997, 101, 6196. (91) Molina, M. J. The probable role of stratospheric ‘ice’ clouds: heterogeneous chemistry of the ‘ozone hole’. In The Chemistry of the Atmosphere: Its Impact on Global Change; Calvert, J. G., Ed.; Blackwell Scientific Publications: Oxford, U.K., 1994. (92) Feynman, R. P. Statistical Mechanics; Benjamin: New York, 1972. (93) Berne, B. J.; Thirumalai, D. Annu. ReV. Phys. Chem. 1986, 37, 401. (94) Tuckerman, M. E. Path Integration via Molecular Dynamics. In Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms; Grotendorst, J., Marx, D., Muramatsu, A., Eds.; John von Neumann Institute for Computing: Ju¨lich, Germany, 2002; Vol. 10, p 269. (95) Schulman, L. S. Techniques and Applications of Path Integration; Wiley & Sons: New York, 1996. (96) Ceperley, D. M. ReV. Mod. Phys. 1995, 67, 279. (97) Chandler, D.; Wolynes, P. G. J. Chem. Phys. 1981, 74, 4078.

5718 J. Phys. Chem. B, Vol. 113, No. 17, 2009 (98) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (99) Hall, W. R.; Berne, B. J. J. Chem. Phys. 1984, 81, 3641. (100) Sprik, M.; Klein, M. L.; Chandler, D. Phys. ReV. B 1985, 31, 4234. (101) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5093. (102) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5106. (103) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6157. (104) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6168. (105) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6184. (106) Jang, S.; Voth, G. A. J. Chem. Phys. 1999, 111, 2357. (107) Jang, S.; Voth, G. A. J. Chem. Phys. 1999, 111, 2371. (108) Reichman, D. R.; Roy, P.-N.; Jang, S.; Voth, G. A. J. Chem. Phys. 2000, 113, 919. (109) Krishna, V.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18953. (110) Jang, S. J. Chem. Phys. 2006, 124, 064107. (111) Paesani, F.; Voth, G. A. J. Chem. Phys. 2008, 129, 194113. (112) Miller, W. H. J. Phys. Chem. A 2001, 105, 2942. (113) Sun, X.; Miller, W. H. J. Chem. Phys. 1999, 110, 6635. (114) Makri, N.; Miller, W. H. J. Chem. Phys. 2002, 116, 9207. (115) Makri, N. J. Phys. Chem. A 2004, 108, 806. (116) Nakayama, A.; Makri, N. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 4230. (117) Liu, J.; Miller, W. H. J. Chem. Phys. 2006, 125, 224104. (118) Liu, J.; Miller, W. H. J. Chem. Phys. 2007, 127, 114506. (119) Liu, J.; Miller, W. H. J. Chem. Phys. 2007, 126, 234110. (120) Liu, J.; Miller, W. H. J. Chem. Phys. 2008, 128, 144511. (121) Liu, J.; Miller, W. H. J. Chem. Phys. 2008, 129, 124111. (122) Reichman, D. R.; Rabani, E. Phys. ReV. Lett. 2001, 87, 265702. (123) Reichman, D. R.; Rabani, E. J. Chem. Phys. 2002, 116, 6271. (124) Poulsen, J. A.; Nyman, G.; Rossky, P. J. J. Chem. Phys. 2003, 119, 12179. (125) Poulsen, J. A.; Nyman, G.; Rossky, P. J. J. Phys. Chem. A 2004, 108, 8743. (126) Poulsen, J. A.; Nyman, G.; Rossky, P. J. J. Phys. Chem. B 2004, 108, 19799. (127) Poulsen, J. A.; Nyman, G.; Rossky, P. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6709. (128) Poulsen, J. A.; Nyman, G.; Rossky, P. J. J. Chem. Theory Comput. 2006, 2, 1482. (129) Poulsen, J. A.; Scheers, J.; Nyman, G.; Rossky, P. J. Phys. ReV. B 2007, 75, 224505. (130) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2004, 121, 3368. (131) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2005, 123, 034102. (132) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2005, 122, 084106. (133) Miller, T. F., III; Manolopoulos, D. E. J. Chem. Phys. 2005, 123, 154504. (134) Miller, T. F., III; Manolopoulos, D. E. J. Chem. Phys. 2005, 122, 184503. (135) Braams, B. J.; Manolopoulos, D. E. J. Chem. Phys. 2006, 125, 124105. (136) Braams, B. J.; Miller, T. F.; Manolopoulos, D. E. Chem. Phys. Lett. 2006, 418, 179. (137) Habershon, S.; Braams, B. J.; Manolopoulos, D. E. J. Chem. Phys. 2007, 127, 174108. (138) Habershon, S.; Fanourgakis, G. S.; Manolopoulos, D. E. J. Chem. Phys. 2008, 129, 074501. (139) Collepardo-Guevara, R.; Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2008, 128, 144502. (140) Horikoshi, A.; Kinugawa, K. J. Chem. Phys. 2003, 119, 4629. (141) Horikoshi, A.; Kinugawa, K. J. Chem. Phys. 2005, 122, 174104. (142) Kuharski, R. A.; Rossky, P. J. J. Chem. Phys. 1985, 82, 5164. (143) Wallqvist, A.; Berne, B. J. Chem. Phys. Lett. 1985, 117, 214. (144) Guillot, B.; Guissani, Y. J. Chem. Phys. 1998, 108, 10162. (145) Hernańdez de la Penã, L.; Kusalik, P. G. J. Chem. Phys. 2004, 121, 5992. (146) Hernańdez de la Penã, L.; Kusalik, P. G. J. Am. Chem. Soc. 2005, 127, 5246. (147) Hernańdez de la Penã, L.; Kusalik, P. G. J. Chem. Phys. 2006, 125, 054512. (148) Miller, T. F., III; Manolopoulos, D. E. J. Chem. Phys. 2005, 123, 154504. (149) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (150) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 120, 300. (151) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515. (152) Lee, H.-S.; Tuckerman, M. E. J. Chem. Phys. 2006, 125, 154507. (153) Lee, H.-S.; Tuckerman, M. E. J. Chem. Phys. 2007, 126, 164501. (154) Chen, B.; Ivanov, I.; Klein, M. L.; Parrinello, M. Phys. ReV. Lett. 2003, 91, 215503. (155) Morrone, J. A.; Car, R. Phys. ReV. Lett. 2008, 101, 017801.

Paesani and Voth (156) Donchev, A. G.; Ozrin, V. D.; Subbotin, M. V.; Tarasov, O. V.; Tarasov, V. I. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 7829. (157) Burnham, C. J.; Li, J. C.; Xantheas, S. S.; Leslie, M. J. Chem. Phys. 1999, 110, 4566. (158) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1500. (159) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115. (160) Fanourgakis, G. S.; Xantheas, S. S. J. Phys. Chem. A 2006, 110, 4100. (161) Stern, H. A.; Berne, B. J. J. Chem. Phys. 2001, 115, 7622. (162) Fanourgakis, G. S.; Schenter, G. K.; Xantheas, S. S. J. Chem. Phys. 2006, 125, 141102. (163) Paesani, F.; Iuchi, S.; Voth, G. A. J. Chem. Phys. 2007, 127, 074506. (164) Burnham, C. J.; Reiter, G. F.; Mayers, J.; Abdul-Redah, T.; Reichert, H.; Dosch, H. Phys. Chem. Chem. Phys. 2006, 8, 3966. (165) Head-Gordon, T.; Johnson, M. E. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 7973. (166) Head-Gordon, T.; Rick, S. W. Phys. Chem. Chem. Phys. 2007, 9, 83. (167) Sorenson, J. M.; Hura, G.; Glaeser, R. M.; Head-Gordon, T. J. Chem. Phys. 2000, 113, 9149. (168) Hura, G.; Russo, D.; Glaeser, R. M.; Head-Gordon, T.; Krack, M.; Parrinello, M. Phys. Chem. Chem. Phys. 2003, 5, 1981. (169) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford: Oxford, U.K., 1987. (170) Yeh, I. C.; Hummer, G. J. Phys. Chem. B 2004, 108, 15873. (171) Hernańdez de la Penã, L.; Razul, M. S. G.; Kusalik, P. G. J. Phys. Chem. A 2005, 109, 7236. (172) Miller, T. F., III; Manolopoulos, D. E. J. Chem. Phys. 2005, 123, 154504. (173) Bertie, J. E.; Lan, Z. Appl. Spectrosc. 1996, 50, 1047. (174) Guillot, B. J. Chem. Phys. 1991, 95, 1543. (175) Corcelli, S. A.; Skinner, J. L. J. Phys. Chem. A 2005, 109, 6154. (176) Lawrence, C. P.; Skinner, J. L. Chem. Phys. Lett. 2003, 369, 472. (177) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 8847. (178) Rey, R.; Moller, K. B.; Hynes, J. T. J. Phys. Chem. A 2002, 106, 11993. (179) Moller, K. B.; Rey, R.; Hynes, J. T. J. Phys. Chem. A 2004, 108, 1275. (180) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698. (181) Eaves, J. D.; Tokmakoff, A.; Geissler, P. L. J. Phys. Chem. A 2005, 109, 9424. (182) Eaves, J. D.; Loparo, J. J.; Fecko, C. J.; Roberts, S. T.; Tokmakoff, A.; Geissler, P. L. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13019. (183) Harder, E.; Eaves, J. D.; Tokmakoff, A.; Berne, B. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 11611. (184) Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2004, 120, 8107. (185) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2004, 121, 8887. (186) Hayashi, T.; Jansen, T. L.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 64. (187) Jansen, T. L.; Hayashi, T.; Zhuang, W.; Mukamel, S. J. Chem. Phys. 2005, 123, 114504. (188) Auer, B.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14215. (189) Moilanen, D. E.; Fenn, E. E.; Lin, Y. S.; Skinner, J. L.; Bagchi, B.; Fayer, M. D. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 5295. (190) Laage, D.; Hynes, J. T. Chem. Phys. Lett. 2006, 433, 80. (191) Tan, H. S.; Piletic, I. R.; Fayer, M. D. J. Chem. Phys. 2005, 122, 174501. (192) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 264. (193) Rezus, Y. L. A.; Bakker, H. J. J. Chem. Phys. 2005, 123, 114502. (194) Winkler, R.; Lindner, J.; Bu¨rsing, H.; Vo¨hringer, P. J. Chem. Phys. 2000, 113, 4674. (195) Fanourgakis, G. S.; Schenter, G. K.; Xantheas, S. S. J. Chem. Phys. 2006, 125. (196) Fenimore, P. W.; Frauenfelder, H.; McMahon, B. H.; Parak, F. G. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 16047. (197) Fenimore, P. W.; Frauenfelder, H.; McMahon, B. H.; Young, R. D. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 14408. (198) Frauenfelder, H.; Fenimore, P. W.; Chen, G.; McMahon, B. H. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 15469. (199) Konermann, L.; Simmons, D. A. Mass Spectrom. ReV. 2003, 22, 1. (200) Schowen, K. B.; Schowen, R. L. Methods Enzymol. 1982, 87, 551. (201) Hummer, G.; Garde, S.; Garcia, A. E.; Pratt, L. R. Chem. Phys. 2000, 258, 349. (202) Case, D. A.; Darden, T. A.; Cheatham, I., T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crowley, M.; Walker, R. C.; Zhang, W.; Merz, K. M.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossva´ry,

Centennial Feature Article I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Yang, L.; Tan, C.; Mongan, J.; Hornak, V.; Cui, G.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; P.A., K. AMBER 10, 10th ed.; University of CaliforniasSan Francisco: San Francisco, CA, 2008. (203) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157. (204) Wang, B.; Merz, K. M. J. Chem. Theory Comput. 2006, 2, 209. (205) Vanicek, J.; Miller, W. H. J. Chem. Phys. 2007, 127, 114309. (206) McIntyre, J. D. E.; Salomon, M. J. Phys. Chem. 1968, 72, 2431. (207) Finlayson-Pitts, B. J.; Pitts, J. N., Jr. Chemistry of the Upper and Lower Atmosphere: Theory, Experiments and Applications; Academic Press: San Diego, CA, 2000. (208) Kroes, G. J. Surf. Sci. 1992, 275, 365. (209) Materer, N.; Starke, U.; Barbieri, A.; VanHove, M. A.; Somorjai, G. A.; Kroes, G. J.; Minot, C. Surf. Sci. 1997, 381, 190. (210) Bolton, K.; Pettersson, J. B. C. J. Phys. Chem. B 2000, 104, 1590. (211) Ikeda-Fukazawa, T.; Kawamura, K. J. Chem. Phys. 2004, 120, 1395.

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5719 (212) Mantz, Y. A.; Geiger, F. M.; Molina, L. T.; Molina, M. J.; Trout, B. L. J. Chem. Phys. 2000, 113, 10733. (213) Conde, M. M.; Vega, C.; Patrykiejew, A. J. Chem. Phys. 2008, 129, 014702. (214) Paesani, F.; Voth, G. A. J. Phys. Chem. C 2008, 112, 324. (215) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: London, 1969. (216) Furukawa, Y.; Yamamoto, M.; Kuroda, T. J. Cryst. Growth 1987, 82, 665. (217) Bluhm, H.; Ogletree, D. F.; Fadley, C. S.; Hussain, Z.; Salmeron, N. J. Phys.: Condens. Matter 2002, 14, L227. (218) Gillen, K. T.; Douglas, D. C.; Hoch, M. J. R. J. Chem. Phys. 1972, 57, 5117. (219) Mills, R. J. Phys. Chem. 1973, 77, 685. (220) Bagchi, B. Chem. ReV. 2005, 105, 3197. (221) Wallqvist, A.; Berne, B. J. J. Chem. Phys. 1993, 97, 13841. (222) Madden, P.; Kivelson, D. AdV. Chem. Phys. 1984, 56, 467. (223) Tan, H. S.; Piletic, I. R.; Fayer, M. D. J. Chem. Phys. 2005, 122.

JP810590C

The Properties of Water: Insights from Quantum Simulations - The

Recommend Documents