Infrared Spectroscopy and Hydrogen-Bond Dynamics of Liquid Water

Sep 1, 2009 - Infrared Spectroscopy and Hydrogen-Bond Dynamics of Liquid Water from Centroid. Molecular Dynamics with an Ab Initio-Based Force Field...
1 downloads 0 Views 5MB Size
13118

J. Phys. Chem. B 2009, 113, 13118–13130

Infrared Spectroscopy and Hydrogen-Bond Dynamics of Liquid Water from Centroid Molecular Dynamics with an Ab Initio-Based Force Field Francesco Paesani,†,‡ Sotiris S. Xantheas,§ 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, and Chemical and Materials Sciences DiVision, Pacific Northwest National Laboratory, 902 Battelle BouleVard, P.O. Box 999, MS K1-83, Richland, Washington 99352 ReceiVed: August 7, 2009

A molecular-level description of the unique properties of hydrogen-bond networks is critical for understanding many fundamental physicochemical processes in aqueous environments. In this Article, a novel simulation approach, combining an ab initio-based force field for water with a quantum treatment of the nuclear motion, is applied to investigate hydrogen-bond dynamics in liquid water with a specific focus on the relationship of these dynamics to vibrational spectroscopy. Linear and nonlinear infrared (IR) spectra are calculated for liquid water, HOD in D2O and HOD in H2O, and discussed in the context of the results obtained using other approaches that have been employed in studies of water dynamics. A comparison between the calculated spectra and the available experimental data yields an overall good agreement, indicating the accuracy of the present simulation approach in describing the properties of liquid water under ambient conditions. Possible improvements on the representation of the underlying water interactions as well as the treatment of the molecular motion at the quantum-mechanical level are also discussed. 1. Introduction The structural and dynamical properties of liquid water play a central role in many processes that occur in Nature. To a significant extent, the unique behavior of water may be attributed to its ability to establish a dynamic hydrogen-bond (H-bond) network in which both the number and strength of the hydrogen bonds continually fluctuate.1 As a consequence, a detailed characterization of its H-bond structure and dynamics is fundamental to the development of a better understanding of the behavior of water in different environments and under different conditions. The generally accepted picture of the local structure of liquid water points to the presence of distorted hydrogen bonds with each H2O molecule being on average 4-fold coordinated in a tetrahedral arrangement. However, a recent interpretation of measured X-ray absorption spectra has suggested that the coordination number in the liquid is actually closer to two with only one donor and one acceptor hydrogen bond for each molecule.2 Subsequent X-ray emission spectra have also been interpreted as emanating from a bimodal distribution whose two peaks have been assigned to tetrahedral and strongly distorted hydrogen-bonded configurations, respectively.3 The proposed structural organization characterized by two distinct motifs with an estimated ratio of highly distorted to tetrahedral-like species of 2:1 is still highly controversial.3-11 Both infrared and Raman spectroscopy are powerful techniques that, by complementing the existing scattering methodologies, can provide direct insights into the structure and dynamics of liquid water. Infrared spectra measured as a * Author to whom correspondence should be addressed. E-mail: [email protected]. † University of Utah. ‡ Present address: Department of Chemistry and Biochemistry, University of CaliforniasSan Diego, 9500 Gilman Drive, La Jolla, CA 92093. § Pacific Northwest National Laboratory.

function of temperature demonstrate that librations (or hindered rotations) play a critical role in the formation of a dynamic H-bond network,1 while the observation of isosbestic points in the Raman spectra suggests the existence of two distinct types of H2O molecules in the liquid phase.12 Similar conclusions were also reached from earlier experiments using ultrafast spectroscopy that revealed the existence of H2O molecules with two different relaxation rates.13 However, subsequent measurements of the temperature dependence of Raman spectra combined with classical molecular dynamics (MD) simulations have shown that the existence of isosbestic points can also be interpreted in terms of a continuous distribution of structures corresponding to increasing distortions around a single-component tetrahedral motif.14,15 Experimentally, the difference in enthalpy between the possible water structures has been estimated to be 2.6 kcal mol-1 (900 cm-1) from the Raman spectra16 and 1.5 kcal mol-1 (500 cm-1) from the X-ray absorption spectra.12 These values, which lie in the range of the energies required for the librational motion, suggest that water molecules in both fundamental and vibrationally excited states may coexist in the liquid phase where they display large amplitude rotational motions. As a consequence, nuclear quantum effects may be expected to play a nonnegligible role in the H-bond dynamics. In this context, the structural properties of liquid H2O and D2O exhibit appreciable differences that are already noticeable under ambient conditions and become more pronounced as the temperature decreases.17-19 Interestingly, resonance Raman spectra probing the structure and thermodynamics of the hydrated electron have demonstrated the preference of the electron for H2O versus D2O in a mixed solvent, which has been explained by the associated change in zero-point energy upon binding.20 Because of intrinsic thermal averaging, linear spectroscopy is not capable of identifying particular signatures that can be directly related to specific groups of water molecules in locally different environments. An inherent problem in the IR spectra

10.1021/jp907648y CCC: $40.75  2009 American Chemical Society Published on Web 09/01/2009

IR Spectroscopy and H-Bond Dynamics of Liquid Water of liquid water is associated with the structure and assignment of the broadband around 3400 cm-1 corresponding to the OH stretching vibrations. There is currently no consensus as to the number or the nature of the underlying absorptions. To this end, it has been previously suggested that the broadband of the OH vibrations in liquid water results from up to five different absorptions that depend on the number of hydrogen bonds,21 their types (single, bifurcated),22 their length,23 and the presence of the bend overtone,24 to name a few. However, the vibrational coupling may also lead to vibrational eigenstates that are delocalized over several molecules, which suggests that the interpretation of the OH band in terms of a single molecule environment is possibly incomplete.25 By contrast, nonlinear infrared spectroscopy overcomes this problem by preparing molecules in spectrally distinct environments and measuring their subsequent relaxation dynamics. Since the OH vibrational frequency of a water molecule is extremely sensitive to the surrounding environment, new insights into the structure and dynamics of the liquid phase have been obtained during the past years from ultrafast infrared spectroscopy.26-72 Narrowband two-color pump-probe measurements of the orientational dynamics of dilute HOD molecules in D2O (which hereafter will be referred to as the HOD:D2O system) show a different frequency behavior from that observed for HOD molecules in H2O (which hereafter will be referred to as the HOD:H2O system).55,59 However, these studies also differed in their experimental conditions. In the case of the HOD:D2O system, the pump pulse was tuned to different frequencies, while, for the HOD:H2O system, the pump pulse was tuned to the center of the absorption band. In order to explain the different relaxation decays, it was initially proposed that tunneling could selectively speed up the reorientation of the most weakly hydrogen-bonded OH groups.59 However, no frequency dependence was observed in IR pump-probe experiments for the HOD:D2O system with 45 fs pulses, which had sufficient bandwidth to span the entire OH absorption band.69 Interestingly, a recent study of the HOD:H2O system with the pump tuned to different frequencies of the absorption band has also found a frequency dependence in the OD orientational dynamics at short times (t < 1.5 ps).43 From a theoretical standpoint, an ideal simulation to probe the H-bond dynamics in liquid water would be the one that treats both electronic and nuclear degrees of freedom at a quantummechanical level. In this context, ab initio molecular dynamics (AIMD) with “on the fly” calculation of the molecular interactions using density functional theory (DFT) has been extensively used to investigate the properties of liquid water.10,73-91 However, it has been shown that the choice of the density functional model is crucial for an accurate description of most of the water properties with different functionals providing different results.88 A recent study has also suggested that some of the commonly used DFT functionals produce a supercooled, overstructured state for water at room temperature, and simulations at elevated temperatures above the corresponding melting points (which were determined at T > 410 K) are needed in order to simulate a liquid phase.92 Furthermore, missing a correct description of the dispersion contributions, most of the functionals used in previous AIMD simulations were in principle not capable of accurately describing the weak intermolecular interactions occurring in the liquid phase. Although several different schemes have been proposed in the literature for improving on the representation of van der Waals interactions within density functional theory (e.g., see refs 93-101), only recently these corrections have been employed in actual simulations demon-

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13119 strating the importance of dispersion contributions for determining both the structural and dynamical properties of the liquid water phase.102 Furthermore, due to the associated computational cost, AIMD simulations with a quantum-mechanical description of the nuclear motion have been limited to the study of structural and thermodynamic properties (e.g., see refs 73 and 83) and no quantum dynamical calculations have been reported so far. On the other hand, valuable insight into the dynamics of dilute HOD molecules in either H2O or D2O have been obtained using different approaches based on classical MD with several empirical water force fields.24,65-67,103-114 These approaches in general provide different results for different thermodynamic and dynamical properties of water.114-116 In order to make a direct comparison with the results obtained from ultrafast spectroscopy experiments, in all of these studies, the OH (OD) stretch Raman and IR lineshapes of HOD were calculated either by using the electric field fluctuations projected along the OH (OD) bond14,15,67 or by employing a two-dimensional map that correlates the OH (OD) vibrational frequencies with the projected electric field. The latter is usually obtained from a combination of classical MD simulations and electronic structure calculations on small water clusters.25,104,109,114,117 Although approximate,9 these approaches have been shown to be very effective in practice. We also note that in previous simulations of the hydrogen-bond dynamics in the HOD:H2O and HOD: D2O systems only the OH (OD) stretching vibration was treated at a quantum-mechanical level, while quantum effects associated with all other degrees of freedom, which may play a nonnegligible role in liquid water,118 were neglected. In this study, we present a detailed investigation of the hydrogen-bond dynamics in liquid water based on a theoretical and computational approach that combines a treatment of the quantum motion of all of the water molecules with an accurate description of the underlying inter- and intramolecular interactions. More specifically, centroid molecular dynamics (CMD)119-126 is employed in conjunction with the recently developed TTM3-F ab initio-based force field for water126 to analyze the dynamics of a single HOD molecule in both light (H2O) and heavy (D2O) water under ambient conditions (T ) 298.15 K). It has been demonstrated that CMD quantitatively accounts for the most important nuclear quantum effects in the condensed phase,127-138 while recent studies have shown that the infrared spectra computed with the TTM3-F model for clusters as well as for bulk water are in good agreement with the available experimental results.126,139 In particular, this study presents a molecular-level analysis of the line shape of the OD (OH) infrared absorption band for an HOD molecule in H2O (D2O) as well as of the corresponding three-pulse vibrational echo signals with an accounting for the nuclear quantum effects on the water structure and dynamics. This in turn allows for a direct and quantitative comparison with the experimental data obtained from ultrafast infrared spectroscopy. The present results are discussed in light of previous classical MD simulations using empirical force fields with a specific focus on the role played by nuclear quantum effects as well as on the description of the underlying water interactions. The paper is organized as follows: In section 2, the theoretical and computational methodologies are presented, while the results are described in section 3. The conclusions are then given in section 4. 2. Theoretical and Computational Methodologies A complete description of the dynamics of an HOD molecule in both liquid H2O and D2O requires a treatment of the molecular

13120

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Paesani et al.

motion at a quantum-mechanical level along with an accurate representation of the underlying Born-Oppenheimer potential energy surface. In order to satisfy both of these requirements, the CMD method119-125 was used in combination with the ab initio-based TTM3-F interaction potential for water.126 Since the CMD method as well as the TTM3-F force field have been previously described in the literature, only a summary of the main concepts and computational details that are specific to the present study are given here. The CMD 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.140 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

(1)

where qc is the imaginary time Feynman path-integral 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

dR1 ... dRPδ(qc - q˜)[

∑ i)1



[

P

fi] exp -β

[

i)1

P



]

∑ P1 V(Ri)

]

(2)

1 dR1 ... dRPδ(qc - q˜) exp -β V(Ri) P i)1

In eq 2, q˜ is the centroid variable corresponding to the center of mass of a ring-polymer with P beads which represents the quantum particle within the path-integral formalism, i.e., P

q˜ )



1 R P i)1 i

(3)

Here, Ri represents the Cartesian position of the ith bead of the ring-polymer and fi is the force acting on it P

fi ) -

∂Φ({Ri}) 1 P i)1 ∂Ri



(4)

In eq 4, Φ({Ri}) corresponds to the potential term of the Hamiltonian in a path-integral representation of the system, i.e., P

Φ({Ri}) )

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

(5)

i)1

In the present study, the centroid variables of each atom were propagated according to eq 1 using the adiabatic separation scheme of ref 122 which results in the so-called adiabatic centroid molecular dynamics (ACMD) method. Briefly, ACMD makes use of a unitary transformation of the Cartesian coordinates of the ring-polymer beads into the corresponding normal modes. Within the normal mode representation, it can easily be shown that the zero-frequency normal mode corresponds to the centroid position of the quantum particle. In order to take

full advantage of the normal-mode transformation, the centroid mass is then set equal to that of the quantum particle, m, while the mass of each nonzero frequency mode is set to µ˜ n ) γ2mωn2, where ωn is the corresponding normal-mode frequency and γ is the so-called adiabaticity parameter.122 The latter has to be chosen to be small enough in order to ensure a sufficiently large separation in time scales between the centroid and all noncentroid motions. It is important to note that, while the centroid variables must strictly follow a Newtonian dynamics, an efficient thermostatting scheme is necessary for all noncentroid modes in order to guarantee a proper sampling of the quantum dispersion associated with each atom. Within ACMD, the centroid thus moves on the potential of mean force generated by the nonzero frequency modes according to the adiabatic theorem.122 The time propagation of the centroid position and momentum provides an efficient means to compute approximate quantum time correlation functions involving linear operators.124,125 In this case, a time correlation function computed within the centroid dynamics formalism represents an approximation to (K) (t) the corresponding Kubo-transformed correlation function CAB (K) (K) ˆ ˆ ˆ ˆ ) 〈A (0)B (t)〉, where A and B are two operators describing the physical observables of interest. If Aˆ is a linear function of (K) (t) is then the position and/or the momentum operators, CAB related to the regular quantum correlation function through the following relationship:

I(ω) )

[ ( ) ]

pβω pβω coth + 1 I(K)(ω) 2 2

(6)

In eq 6, I(K)(ω) and I(ω) are the Fourier transforms of the Kubotransformed and regular quantum time correlation functions, respectively. By contrast, in the case of nonlinear operators, the situation is more complicated, which has led to the development of different numerical schemes.141-144 In all ACMD simulations discussed in the following section, the interactions between water molecules were described using the ab initio-based TTM3-F potential,126 which corresponds to the latest refinement of the TTM family of force fields for water developed in recent years by Xantheas and co-workers.126,145-149 The family of TTM models, which is based on the Thole-type approach,150 was parametrized using the results of high-level ab initio calculations of the water dimer potential energy surface145 and validated from the energetics of similar results for small water clusters.146,151-153 Within this framework, nonadditive contributions are included through an induction scheme in which the atomic dipoles are determined selfconsistently. This specific parametrization approach was chosen because it was demonstrated that for water clusters the twobody pairwise additive interaction corresponds to more than 80% of the total energy154 and the three-body nonadditive term is the next most important contribution, with higher-order terms being almost negligible.155,156 In the flexible TTM models, the intramolecular vibrations of the water molecules are described using the gas phase monomer potential energy surface developed by Partridge and Schwenke157 in combination with a nonlinear dipole moment surface (DMS). The latter was shown to be a critical component for reproducing the measured increase of the H-O-H angle in water clusters, liquid water, and ice from its gas phase monomer value.158 In the latest version of the model (TTM3-F), the gas phase monomer DMS was modified in order to account for the correct dissociation limits in the condensed phase and was shown to reproduce the experimental

IR Spectroscopy and H-Bond Dynamics of Liquid Water infrared spectra of both water clusters and liquid water with a high degree of accuracy.126 Since the TTM parametrization provides an approximate but accurate representation of the actual multidimensional BornOppenheimer potential energy surface, it is more appropriate to use this model in conjunction with quantum rather than classical simulation schemes. In this regard, path-integral molecular dynamics and CMD calculations of the properties of both liquid water and ice carried out with a previous version of the TTM model (TTM2.1-F) showed a significantly better agreement with the experimental results than their classical counterparts.130,159 Furthermore, since polarization effects are explicitly included in the TTM parametrization, the atomic charges of each water molecule vary instantaneously in response to any structural modification of the local environment. As a consequence, the TTM model does not assume any a priori structure of the liquid but is instead capable, in principle, of capturing the occurrence of possible asymmetries in the H-bond network. In this study, the ACMD simulations were carried out at a temperature of T ) 298.15 K for a system consisting of a single HOD molecule in a cubic box with other 127 H2O or D2O molecules at the experimental density of light or heavy water, respectively. The short-range interactions were truncated at an atom-atom distance of 7.8 Å, while the electrostatic interactions were treated using the Ewald sum. For both the HOD:H2O and HOD:D2O systems, the results reported in the following section were obtained by averaging over 25 ACMD independent trajectories for a total simulation time of ∼2 ns. All calculations were performed with P ) 32 beads, and the adiabaticity parameter was set to γ ) 0.4, which was found to be small enough to guarantee a sufficiently large separation in time scales between the centroid and noncentroid motions. The ACMD equations of motion were propagated according to the velocityVerlet algorithm with a time step of ∆t ) 0.02 fs. In order to provide a direct comparison with the available experimental data from pump-probe infrared measurements, the instantaneous OD (OH) frequency corresponding to the vibrational transition between the ground and first excited state (ω10) and between the first two excited states (ω21) of the HOD molecule in H2O (D2O) was also computed during each ACMD trajectory from a numerical solution of the Schro¨dinger equation for the corresponding one-dimensional OD (OH) oscillator using Numerov’s method.160 Following the procedure described in ref 108 for any given configuration of the HOD:H2O (HOD:H2O) system, the OD (OH) Born-Oppenheimer potential energy curve was computed by stretching the OD (OH) bond while keeping the positions of all other atoms fixed. The resulting potential curve was then used to compute the corresponding vibrational frequency as well as the vibrational wave function of the ground and first two excited states. The vibrational transition moments were obtained within the same computational scheme from the numerical integration over the product of the molecular dipole moment and the appropriate vibrational wave functions. In all calculations, the intramolecular motion of the HOD molecule was represented in terms of the molecular Hamiltonian described by the TTM3-F model. Therefore, in contrast to all previous simulations of the ultrafast hydrogenbond dynamics, both inter- and intramolecular interactions of the water molecules were consistently represented on the same footing by the TTM3-F force field. 3. Results and Discussion 3.1. Linear Spectrum of Bulk H2O. In order to assess the accuracy of the approach discussed in the previous section in

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13121

Figure 1. Infrared spectrum (red line) of light (bottom panel) and heavy (top panel) water at T ) 298.15 K calculated from eq 7 using centroid molecular dynamics. Also shown as black lines are the experimental data from refs 162 and 163, respectively.

reproducing the spectroscopic properties of liquid water and their relationship to the underlying hydrogen-bond dynamics, the infrared absorption coefficient for both light and heavy water at T ) 298.15 K was computed from the expression R(ω) )

πω [ 3Vpcn(ω) ](1 - e

-βpω

)

1 2π



+∞

-∞

dt e-iωt〈µ(0) · µ(t)〉

(7) where n(ω) is the refractive index of the system at frequency ω, c is the speed of light, and 〈µ(0) · µ(t)〉 is the quantum dipole autocorrelation function, which, in principle, is a (weakly) nonlinear function of the atom positions. However, although a direct calculation of the infrared absorption coefficient using the method proposed in ref 144, which combines ACMD with a maximum entropy analytic continuation (MEAC) of the relevant imaginary-time correlation function, proved difficult due to the slow convergence of the imaginary-time dipole autocorrelation function as well as to the complexity of the spectral function itself, a careful inspection of the (strongly) nonlinear HOD orientational correlation function showed that the differences between the ACMD and ACMD + MEAC results are in this case negligible. This is in agreement with the results of ref 144, showing that, in the high-temperature regime (that is the case of the present simulations), the ACMD correlation functions with a classical description of the quantum operators represent very good approximations to the actual nonlinear quantum correlation functions even for highly anharmonic systems.144 It is also important to note that CMD was previously shown to reproduce rather accurately the vibrational red shift in hydrogen-bonded systems.161 Therefore, in the present study, the quantum dipole autocorrelation function was approximated by the corresponding ACMD autocorrelation function. Figure 1 shows a comparison between the present calculations of the infrared spectrum of H2O (bottom panel) and D2O (top panel) along with the corresponding available experimental data.162,163 Both spectra display similar features with the

13122

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Figure 2. Normalized infrared lineshapes for the OD (bottom panel) and OH (top panel) stretching vibration of an HOD molecule in H2O and D2O, respectively. Red line, present results obtained from eq 8 using centroid molecular dynamics; black line, experimental data from refs 162 and 163, respectively.

absorption frequencies for heavy water being significantly redshifted as a consequence of the smaller zero-point energy associated with the motion of D2O molecules. Overall, good agreement is also found with the corresponding experimental results. In particular, the calculated intensity and line shape of the absorption bands associated with the bending vibrations (1550-1600 cm-1 for H2O and 1200-1250 cm-1 for D2O) as well as with the stretching vibrations (3000-3700 cm-1 for H2O and 2000-2800 cm-1 for D2O) are reproduced, although their positions appear to be slightly shifted (∼100 cm-1) toward lower frequencies. Interestingly, the corresponding classical infrared spectrum calculated for liquid water shows instead an apparently better agreement with the experimental results in the OH stretching region.126 Somewhat less satisfactory is the agreement with the available experimental data at frequencies below 1000 cm-1, which correspond to the absorption band associated with the librational motion of the water molecules. In this case, the simulated spectrum predicts a narrower band that implies the absence of relatively high-frequency hindered rotations. The calculated spectra also exhibit the shoulder at ∼200 cm-1, which arises from a dipole-induced dipole mechanism,164 albeit its shape is not well reproduced. A similar infrared spectrum for liquid H2O has recently been reported using the PA-CMD method.139 3.2. Infrared Lineshapes for HOD in H2O and D2O. As mentioned in the Introduction, the OH (OD) vibrational frequency of a HOD molecule in D2O (H2O) is extremely sensitive to the surrounding environment. Therefore, a molecular-level insight into the structural and dynamical properties of liquid water can be obtained from a detailed analysis of the dependence of this frequency on the local structure of the hydrogen-bond network. In order to make a direct comparison with the results from experiments of ultrafast infrared spectroscopy, the OD (bottom panel) and OH (top panel) lineshapes, I(ω), calculated for an HOD molecule in both H2O and D2O are shown in Figure 2

Paesani et al.

Figure 3. Comparison between the calculated static frequency distribution (red line) and the infrared line shape obtained from eq 8 using centroid molecular dynamics (black line). Bottom panel, results for the HOD:H2O system; top panel, results for the HOD:D2O systems.

along with the experimental results.29,68 The theoretical lineshapes were obtained from the corresponding ACMD dipole autocorrelation functions with a classical representation of the dipole operator according to the following expression

IACMD(ω) ∼

∫-∞+∞ dt e-iωt〈µ(0) · µ(t)〉

(8)

As expected from the analysis of the infrared spectra of light and heavy water (Figure 1), the center of both OH and OD lines is shifted by ∼100 cm-1 with respect to the experimental values. However, their shapes are found to be in fairly good agreement, although a slightly too long tail appears at higher frequencies which is apparently at odds with the experimental results. The calculated full width at half-maximum (fwhm) is ∼140 cm-1 for the OD line and ∼220 cm-1 for the OH line which compare quite well with the corresponding experimental values of ∼160 and ∼250 cm-1. In principle, the infrared lineshapes contain all the information regarding the structure and dynamics of liquid water. However, because of the intrinsic thermal averaging, this information is hidden underneath the relatively broad spectral envelope, which makes difficult, if not impossible, the identification of particular signatures specific to well-defined local environments. In order to overcome this difficulty and to make a direct connection between infrared spectra and H-bond configurations, the distributions of both OD and OH instantaneous frequencies corresponding to transitions between the ground and first vibrational states of the HOD molecule were calculated according to the scheme discussed in section 2 using the ACMD trajectories for the HOD:H2O and HOD:D2O systems, respectively. The two frequency distributions P(ω) are shown in the bottom (for OD) and top (for OH) panels of Figure 3. Both P(ω) were obtained by binning ∼4 × 106 frequencies representing as many independent configurations extracted from the corresponding ACMD trajectories. The calculated spectral lines

IR Spectroscopy and H-Bond Dynamics of Liquid Water

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13123

and frequency distributions display similar shapes, although the latter are significantly broader which may likely be due to the neglect of motional narrowing and possible non-Condon effects in the calculation of the static frequency distributions. In order to assess the relative importance of such effects, the OD and OH spectral lines were also computed using mixed quantumclassical approaches within various degrees of approximation. As described in detail in ref 165, the OD (OH) spectral line for the 1-0 vibrational transition of dilute HOD in H2O (D2O) can be expressed in the classical limit as

∫-∞+∞ dt e-iωt〈µ10(0) · µ10(t)e[i ∫ dτ ω t

Incon(ω) ∼

〉e-t/2T1 (9)

10(τ)]

0

where µ10(t) is the vibrational transition moment associated to the transition between the ground and excited state of the OH (OD) bond, ω10(t) is the instantaneous vibrational frequency at time t, T1 is the lifetime of the first excited state (T1 ) 1450 fs for the OD vibration of HOD in H2O, and T1 ) 700 fs for the OH vibration of HOD in D2O114), and the brackets denote a classical equilibrium average. Within the Condon approximation, it is then assumed that the magnitude of the vibrational transition moment is independent of the surrounding environment. Considering that the orientational dynamics of a water molecule is significantly slower than its stretching vibrational motion, the expression for the spectral line in eq 9 becomes

Icon(ω) ∼

∫-∞

+∞

dt e-iωt〈e[i

∫ dτ ω t

10(τ)]

0

〉e-t/2T1

Figure 4. Frequency time correlation function, eq 13, for the OD and OH 1-0 vibrational transitions of an HOD molecule in H2O (bottom panel) and D2O (top panel), respectively. Red lines, present results; black lines, experimental data from refs 64, 67, and 171 (a, b, and c, respectively) and from ref 29 (d).

(10)

Equation 10 can be further simplified by performing a cumulant expansion truncated at second order, which is exact in the limit of a Gaussian distribution of the frequency fluctuations, leading to the following expression

Icum(ω) ∼

∫-∞+∞ dt e-i(ω-〈ω

10〉)t

e-g10,10(t)e-t/2T1

(11)

Here, 〈ω10〉 is the average frequency and

gij,kl(t) )

∫0t dτ(t - τ)cω ω (τ) ij kl

(12)

cωijωkl(t) with i ) 1, j ) 0, k ) 1, and l ) 0 is the frequency time correlation function (FTCF) involving the ground and first excited vibrational states

cωijωkl(t) ) 〈δωij(t)δωkl(0)〉

(13)

with δωij(t) ) ωij(t) - 〈ωij〉. A comparison between the calculated and several experimentally derived frequency time correlation functions is shown in Figure 4 for HOD in D2O (top panel) and HOD in H2O (bottom panel). In the latter case, the calculated FTCF appears to have a faster decay than its experimental counterpart. However, it is difficult to give any conclusive statement about the overall accuracy of the present results considering the clear differences among the available experimental data for the HOD:D2O system which have recently been shown to reflect qualitatively different broadening mechanisms ranging from modestly inhomogeneous to fully homogeneous ones.166

Figure 5. Comparison between several mixed quantum-classical approximations to the infrared lineshape and the corresponding static frequency distribution (black dotted line) and the infrared line shape obtained from eq 8 using centroid molecular dynamics (black dashed line). Red line, cumulant approximation from eq 11; green line, Condon approximation from eq 10; non-Condon approximation from eq 9. Bottom panel, results for the HOD:H2O system; top panel, results for the HOD:D2O systems.

The IR lineshapes for all three mixed quantum-classical approximations are shown in Figure 5 for the OH (top panel) and OD stretch (bottom panel) of HOD in D2O and H2O, respectively. Also shown are the corresponding frequency

13124

J. Phys. Chem. B, Vol. 113, No. 39, 2009

distributions and ACMD lineshapes. For both the HOD:D2O and HOD:H2O systems, all mixed quantum-classical lineshapes are appreciably less broad than the corresponding frequency distributions, a fact which suggests that motional narrowing effects in liquid water are indeed important. The comparison in Figure 5 also indicates that, for bulk water configurations obtained from ACMD simulations with the TTM3-F force field, the Condon and non-Condon approximations result in spectral lines that are slightly red-shifted and narrower by ∼15-20 cm-1 than the corresponding lines calculated within the cumulant approximation. This is in qualitative agreement with the results reported in ref 167 for the SPC/E model. However, the overall shape of the infrared line changes only marginally going from the cumulant to Condon and non-Condon approximations with the last two being, nevertheless, capable of better reproducing the small asymmetry found on the blue side of the ACMD line shape. These findings are not completely unexpected, since, for Gaussian-like frequency distributions such as those shown in Figure 2, the cumulant approximation already provides a very accurate description of the spectral lines. A comparison with the analysis of ref 167 however suggests that non-Condon effects may be strongly dependent on the water model employed in the simulations. Indeed, it was shown in ref 167 that for water configurations obtained with the SPC/E model (which, contrary to the TTM3-F model, provides significantly more asymmetric frequency distributions) the non-Condon effects play a major role in an accurate description of the actual infrared line shape. Figure 5 also illustrates the remarkably good agreement between the ACMD and all three mixed quantum-classical lineshapes. Such a level of agreement clearly demonstrates that the blue shift of about 100 cm-1 obtained for the OH stretching band using centroid molecular dynamics with respect to the results of classical MD simulations139 is a direct result of the quantization of the nuclear vibrations as well as of the specific features of the underlying Born-Oppenheimer surface and not an artifact of the ACMD method. As mentioned in the Introduction, measurements of Raman and X-ray absorption spectra indicate that water molecules in both ground and excited librational states may coexist in the liquid phase at room temperature, which implies that nuclear quantum effects may play a non-negligible role in the water dynamics. Since the ACMD method is capable of accurately describing the librational motion at a quantum-mechanical level, the present analysis provides a unique opportunity to quantitatively assess the importance of nuclear quantum effects on the dynamics of the H-bond network. For this purpose, additional classical MD simulations were performed for the HOD:H2O system with the TTM3-F model. A comparison between the OD frequency distributions obtained from ACMD and classical trajectories is shown in Figure 6. The two distributions clearly display similar shapes. However, the ACMD results are redshifted by about 10-15 cm-1 as a consequence of the explicit quantization of the motion of all water molecules. This comparison thus indicates that nuclear quantum effects are important for a quantitative description of the water dynamics, especially when the underlying molecular interactions are described from a first-principles-based force field as in the present study. However, the similarity between the ACMD and classical frequency distributions in Figure 6 also implies that, under ambient conditions, such effects do not qualitatively modify the behavior of the H-bond network, which is in agreement with a previous analysis based on the temperature dependence of diffusive and relaxation processes in bulk water.130

Paesani et al.

Figure 6. Comparison between the frequency distributions computed for the HOD:H2O system using the ACMD method (red line) and classical MD simulations (black line).

3.3. Relationship between the Electric Field and Infrared Lineshapes. A molecular-level description of the H-bond dynamics in liquid water has recently started to emerge from the analysis of the electric field fluctuations projected along the direction of the OH (OD) bond of the HOD molecule in D2O (H2O) which can easily be calculated using molecular dynamics simulations. In particular, it was shown that the OH (OD) vibrational frequency of the HOD molecule is correlated with the electric field parallel to the OH (OD) bond, which is generated by the surrounding molecules and acts on the H (D) atom.168 Different computational schemes have been developed in the past few years exploiting this correlation in order to calculate infrared and Raman spectra that have then been interpreted in terms of the corresponding molecular configurations obtained from classical MD simulations.14,25,65,67,104,109,114,117,169,170 The electric field fluctuations obtained from classical molecular dynamics simulations with the SPC/E water force field have been directly used to analyze two-dimensional infrared spectra as well as to study the Raman spectra of liquid water.14,65,67,170 Two-dimensional maps that correlate the electric field projected along the OH or OD bonds of HOD with the corresponding vibrational frequencies and dipole moments have also been constructed on the basis of DFT calculations for small water clusters extracted from classical molecular dynamics trajectories. Due to the computational cost associated with the electronic structure calculations, only a limited number of water clusters (∼100-1000) has been considered. This approach has been employed with success in the study of linear and nonlinear infrared spectra as well as of Raman spectra of liquid water.25,104,109,114,117,169 In order to provide more detailed and quantitative insights, the relationship between the electric field and the spectroscopic signatures was investigated using configurations from the ACMD trajectories with the TTM3-F force field for both the HOD:D2O and HOD:H2O systems. The current approach thus provides the unique opportunity to directly compare the electric field distributions not only with the corresponding vibrational frequency distributions but also with the (quantum) spectral lines calculated within the CMD formalism. Furthermore, since both molecular configurations and OH (OD) vibrational frequencies are obtained consistently employing the TTM3-F water model without relying on additional ab initio calculations, this affords the possibility to consider a large number of water configurations

IR Spectroscopy and H-Bond Dynamics of Liquid Water

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13125

Figure 7. Two-dimensional maps correlating the electric field along the OD (OH) bond of the HOD molecule and the corresponding OD (OH) vibrational frequency. Bottom panel, results for the HOD:H2O system; top panel, results for the HOD:D2O system. White dotted line: linear regression of the data.

for the entire system, which significantly improves the statistical representation of possible water structures as well as provides a more accurate description of the molecular interactions in the bulk phase. For the purpose of the present analysis, the component of the electric field (E-field) along the OD (OH) bond of the HOD molecule due to the surrounding H2O (D2O) molecules was computed along the ACMD trajectories with the TTM3-F water model. The electric field was then directly correlated to the instantaneous OD (OH) vibrational frequencies calculated for the same configuration as described in section 2, which allowed us to construct a two-dimensional map from ∼4 × 106 independent configurations. The resulting maps, which correlate the electric field with the OD and OH frequencies for both the HOD:H2O and HOD:D2O systems, are shown in Figure 7. Also indicated as dotted white lines are the best linear interpolations computed for both sets of data. Clearly, the present data show a strong relationship between the electric field and the OH (OD) vibrational frequency with a correlation coefficient (assuming a linear correlation) of approximately -0.9 for both the HOD: H2O and HOD:D2O systems, which is similar to the values reported in the literature for the SPC/E water model.109 An analysis of a significantly smaller number of randomly chosen water configurations (i.e., 100-1000, which is comparable to the ensembles used in previous classical MD/electronic structure calculations)109,117 provides an essentially identical fit, which again demonstrates the intrinsic relationship between the OH (OD) vibrational frequency and the component of the electric field along the corresponding bond. The two-dimensional maps of Figure 7 also indicate that the range of the calculated E-field values is very similar to those obtained from previous DFT calculations, a result that reinforces the notion that the TTM3-F water model provides an accurate representation of the underlying multidimensional BornOppenheimer potential energy surface for liquid water. However, a more detailed inspection reveals that the TTM3-F model predicts a somewhat shifted and narrower range of E-field values lacking contributions below 0.02 au. Unfortunately, at this stage, it is difficult to provide a quantitative assessment of this difference that may likely be due to one or a combination of several different factors, such as small inaccuracies in the parametrization of TTM3-F water force field and intrinsic

Figure 8. Comparison between the electric field distribution (red line) and the corresponding static frequency distribution (black dotted line) and infrared line shape obtained from eq 8 using centroid molecular dynamics (black dashed line). Bottom panel, results for the HOD:H2O system; top panel, results for the HOD:D2O systems.

deficiencies of the DFT functionals. Therefore, a more extensive investigation of the electric field in liquid water at a higher ab initio level of theory is necessary in order to address this issue and will also play a central role in future refinements of the TTM3-F force field. The linear fits shown in the two-dimensional maps were then used to directly investigate the relationship between the electric field fluctuations and the infrared spectral lines. This comparison is reported in Figure 8 which shows the E-field distributions obtained from the linear regression of the data in Figure 7 along with the corresponding frequency distributions and ACMD spectral lines for both the HOD:D2O (top panel) and HOD: H2O (bottom panel) systems. Evidently, all sets of calculations are in qualitative agreement with each other. As expected from the high correlation of the data in the two-dimensional maps of Figure 7, an essentially perfect overlap is found between the electric field and frequency distributions, with both being slightly broader than the ACMD line shape. These results thus clearly demonstrate that the E-field fluctuations provide a fairly reasonable approximation to the OH (OD) infrared spectral line, although inclusion of dynamical effects (see section 3.2) appears to be essential in order to reproduce the correct lineshape. In order to investigate the sensitivity of the OH (OD) vibrational frequency on the surrounding environment, twodimensional maps of E-field vs frequency were also computed for water clusters extracted from the ACMD trajectories containing an increasing number of H2O molecules (N ) 4, 10, and 30) around the single HOD molecule. The results shown in Figure 9 clearly demonstrate that the strong correlation between the vibrational frequency and the electric field is already apparent for small water clusters. Furthermore, the resulting twodimensional maps display a high level of similarity with that obtained from bulk simulations (Figure 6). Assuming a perfect linear correlation (the correlation coefficient in all three cases is >0.93), the E-field fluctuations were then used to compute

13126

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Paesani et al. from ultrafast infrared spectroscopy, the three-pulse echo signal was also analyzed. Experimentally, this signal is generated from a sequence of three subsequent laser pulses with specific wave vectors k1, k2, and k3 and time delays t1 (between the first two pulses) and t2 (between the second two pulses). The resulting signal is then detected at a time t3 after the third pulse in the background-free directions kr ) -k1 + k2 + k3 and knr ) k1 - k2 + k3, which correspond to the rephasing and nonrephasing signals, respectively. Since the calculation of multitime correlation functions within the CMD formalism is still problematic, only the mixed quantum-classical approach was employed in the present study. Within this framework, the three-pulse echo spectrum is obtained from the double Fourier-Laplace transforms over t1 and t3 of the rephrasing, Rr(t3, t2, t1), and nonrephasing, Rnr(t3, t2, t1), thirdorder response functions

Figure 9. Two-dimensional maps correlating the electric field along the OD bond of the HOD molecule and the corresponding OD vibrational frequency computed for HOD(H2O)N clusters with N ) 4, 10, and 30.

I(ω1, t2, ω3) ≡ Re

∫0∞ dt1 ∫0∞ dt3[ei(-ω t +ω t )Rr(t3, t2, t1) + 11

33

ei(ω1t1+ω3t3)Rnr(t3, t2, t1)] (14) For the purposes of the present analysis and also considering that non-Condon effects have been shown in section 3.2 to play a marginal role in determining the infrared lineshapes with the TTM3-F force field, the three-pulse echo signal was only computed within the cumulant approximation. A more detailed analysis of non-Condon effects on the nonlinear spectroscopy of water modeled with the TTM3-F force field will be the subject of future work. The nonlinear response functions in eq 14 were then calculated from the corresponding cumulant expressions as

Rr(t3, t2, t1) ) 2|〈µ10〉| 4 exp[i〈ω10〉(t1 - t3) - g10,10(t1) + g10,10(t2) - g10,10(t3) - g10,10(t1 + t2) - g10,10(t1 + t3) + g10,10(t1 + t2 + t3)] - |〈µ10〉| 2 |〈µ21〉| 2 exp[i〈ω10〉t1 i〈ω21〉t3 - g10,10(t1) + g10,21(t2) - g21,21(t3) g10,21(t1 + t2) - g10,21(t1 + t3) + g10,21(t1 + t2 + t3)] (15)

Figure 10. Comparison of the electric field distributions computed for HOD(H2O)N clusters with N ) 4, 10, and 30 with that obtained in bulk simulations.

the corresponding OD infrared lineshapes which are compared in Figure 10 with the results obtained from the bulk simulations. While appreciable differences are present between the bulk line shape and that obtained from the HOD(H2O)4 clusters, from the comparison in Figure 10, it is evident that 10 water molecules are already sufficient to reproduce the most important features of the bulk line shape. A slightly better agreement is obtained with clusters containing 30 water molecules (corresponding roughly to the full second solvation shell) around HOD. Similar agreement was also found in a direct comparison of the corresponding frequency distributions (not shown), which is not surprising given the almost perfect correlation obtained in the present analysis of the ACMD trajectories with the TTM3-F model between the vibrational frequency and electric field. These results thus demonstrate that the OH (OD) vibrational frequency of HOD is barely affected by molecules beyond the second solvation shell. Furthermore, the great similarity between the bulk and cluster results also sheds new light into the physical ground behind the success of twodimensional maps in reproducing the infrared lineshapes. 3.4. Two-Dimensional Infrared Spectra. In order to make a more direct comparison with the experimental data available

and

Rnr(t3, t2, t1) ) 2|〈µ10〉| 4 exp[-i〈ω10〉(t1 + t3) - g10,10(t1) g10,10(t2) - g10,10(t3) + g10,10(t1 + t2) + g10,10(t1 + t3) g10,10(t1 + t2 + t3)] - |〈µ10〉| 2 |〈µ21〉| 2 exp[-i〈ω10〉t1 i〈ω21〉t3 - g10,10(t1) - g10,21(t2) - g21,21(t3) + g10,21(t1 + t2) + g10,21(t1 + t3) - g10,21(t1 + t2 + t3)] (16) Here, 〈µij〉 is the average transition moment between vibrational states i and j, while gil,kl(t) is given by eq 12 with the frequency time correlation function involving the vibrational states i, j, k, and l. The relevant 〈µij〉 and gil,kl(t) were computed by solving the one-dimensional Schro¨dinger equation for the OD (OH) stretching motion of the HOD molecule in H2O (D2O) using the same ACMD molecular configurations as in the case of the linear infrared lineshapes in section 3.2. A comparison between the three-pulse vibrational echo spectra for the HOD:D2O (top panel) and HOD:H2O (bottom panel) systems calculated at different values for the time delay

IR Spectroscopy and H-Bond Dynamics of Liquid Water

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13127

Figure 11. Two-pulse vibrational echo spectra computed for several delay times t2 within the cumulant approximation, eqs 14-16. Bottom panel, results for the HOD:H2O system; top panel, results for the HOD:D2O systems.

t2 is shown in Figure 11, where the positive values (warm colors) correspond to the 1-0 resonance while negative values (cold colors) correspond to the 2-1 resonance. Clearly, the two sets of results display very similar time dependences, indicating that the hydrogen-bond dynamics is similar in either light or heavy water under ambient conditions. This similarity is even more evident in Figure 12, which shows a comparison between the calculated and experimental nodal slope. The latter is the slope of a straight line separating the positive region of the 1-0 resonance from the negative region of the 2-1 resonance in the three-pulse vibrational echo spectra of Figure 11. For the HOD:D2O system, the present results are in reasonably good agreement with the available experimental data especially at short times, although they apparently miss the oscillation at about 150 fs. By contrast, for the HOD:H2O system, the calculated nodal slope is appreciably smaller than the experimental counterpart and also displays a different decay behavior. A comparison with the results reported in ref 114 for the nodal slope indicates that the present calculations with the TTM3-F force field show a similar agreement with the experimental data as the SPC/E model. In order to better assess the reliability of the TTM3-F force field in describing the dynamics in liquid water, the dynamic line width (DLW) was also computed for the HOD:H2O system. The DLW is defined as the fwhm of a slice of the 1-0 resonance in the three-pulse vibrational echo spectrum taken through its maximum. A comparison between the values extracted from the two-dimensional IR spectra in Figure 11 and the experimental data from ref 28 is shown in Figure 13. In the ideal case of infinitely short pulses (impulsive limit), the DLW approaches the fwhm of the linear spectrum as the delay time t2 increases. However, in the real case of finite pulses, the DLW approaches a slightly narrower value that is dependent on the pulse profiles. Therefore, since the calculated infrared spectrum

Figure 12. Comparison between the calculated (red) and experimental (black) nodal slope. The experimental data for the HOD:D2O system are from ref 65, while those for the HOD:D2O system are from ref 28.

is slightly narrower than its experimental counterpart (Figure 2) and, moreover, the cumulant approximation also underestimates the line width of the ACMD spectral line (Figure 5), the DLW calculated from the two-pulse vibrational echo spectra appears to be consistently smaller than the corresponding experimental values. However, it is important to note that the

13128

J. Phys. Chem. B, Vol. 113, No. 39, 2009

Figure 13. Comparison between the calculated (red) and experimental (black) dynamic line width for the HOD:H2O system. The experimental data are from ref 28.

time dependence as a function of t2 is well reproduced, which reinforces the notion that the TTM3-F force field in combination with a quantum treatment of the nuclear motion within the CMD formalism is capable of describing reasonably well the H-bond dynamics of liquid water under ambient conditions. Importantly, since the water configurations used in this study to calculate the three-pulse vibrational echo signals were obtained from ACMD trajectories, which explicitly take into account nuclear quantum effects, the similarity between the twodimensional spectra for the HOD:D2O and HOD:H2O systems implies that at ambient temperature tunneling plays a negligible role in the hydrogen-bond dynamics as inferred from the frequency distributions computed in section 3.2 at a quantum and classical level of theory, and in agreement with previous analyses of the water orientational dynamics.130 4. Conclusions In this work, a molecular-level investigation of the IR spectroscopy of liquid water under ambient conditions has been carried out employing a simulation approach that combines CMD for treating the nuclear motion at a quantum-mechanical level with an accurate representation of the underlying multidimensional Born-Oppenheimer surface described by the ab initio-based TTM3-F water interaction potential. The linear spectra computed for both light and heavy water from the corresponding (quantum) dipole autocorrelation functions are found in good agreement with the experimental results. The calculated bending and stretching absorption bands are shifted by ∼100 cm-1 toward lower frequencies, while the librational band below 1000 cm-1 appears to be somewhat too narrow. In order to make a direct connection with the experimental ultrafast spectroscopic measurements, the OH and OD frequency distributions and infrared lineshapes for an HOD molecule in D2O and H2O, respectively, have also been analyzed using different theoretical approaches. The frequency distributions have been found to be significantly broader than the corresponding lineshapes obtained from the dipole autocorrelation functions within the CMD framework, indicating that motional narrowing effects are important for the infrared spectroscopy in liquid water. This has been confirmed by a detailed analysis of several mixed quantum-classical approximations to the vibrational lineshapes. An excellent agreement is found between the ACMD and mixed quantum-classical spectral lines, a fact

Paesani et al. that clearly demonstrates that centroid molecular dynamics represents an accurate and efficient method for a direct calculation of infrared spectra in the condensed phase from dipole autocorrelation functions. The relationship between the magnitude of the electric field acting on the OH (OD) bond of the HOD molecule and the corresponding vibrational frequencies has also been investigated considering ∼4 × 106 independent (quantum) molecular configurations for both the HOD:D2O and HOD:H2O systems. Twodimensional maps that relate the electric field value to the vibrational frequency for a specific configuration show an intrinsic strong correlation. Assuming a perfect linear correlation, the OH (OD) infrared line shape has been thus also calculated from the corresponding distribution of the electric field fluctuations. The resulting line shape is found in excellent agreement with the calculated frequency distributions, suggesting that the electric field distributions provide a reasonable approximation to the actual OH (OD) infrared spectral line, although dynamical effects are necessary to recover the correct lineshape. A comparison with the DFT-based results reported in the literature indicates that the range of electric field values predicted by the TTM3 force field is slightly narrower and shifted toward higher values, which is very likely responsible for the red shift found in the infrared spectrum calculated using ACMD. A detailed analysis of the effects of the surrounding environment demonstrates that the OH (OD) vibrational frequency of HOD is barely affected by water molecules beyond the second solvation shell, which provides a new insight into the physical basis behind the success of two-dimensional maps in reproducing the infrared lineshapes. Finally, the analysis of two-pulse vibrational echo spectra has reinforced the notion that quantum simulations within the CMD formalism with the ab initio-based TTM3-F force field provide a reasonably accurate description of the H-bond dynamics in liquid water under ambient conditions. Future work will focus on the refinement of the TTM3-F force field in order to obtain better agreement with the experimental results for linear and nonlinear infrared spectra, as well as on the development of a theoretical approach for computing nonlinear response functions within the centroid dynamics framework. Acknowledgment. This research was supported by the Office of Naval Research through grant N00014-05-1-0457 and by the Division of Chemical Sciences, Biosciences and Geosciences, US Department of Energy. Battelle operates the Pacific Northwest National Laboratory for the US Department of Energy. The authors thank Prof. Jim Skinner and Dr. Piotr Pieniazek for valuable discussions on the hydrogen-bond dynamics in liquid water and for sharing their data for the SPC/E force field. References and Notes (1) Mare´chal, Y. The hydrogen bond and the water molecule: The physics and chemistry of water, aqueous and bio-media; Elsevier: Amsterdam, The Netherlands, 2006. (2) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; Na¨slund, L. Å.; Hirsh, T. K.; Ojama¨e, L.; Glatzel, P.; Pettersson, L. G. M.; Nilsson, A. Science 2004, 304, 995–999. (3) Tokushima, T.; Harada, Y.; Takahashi, O.; Senba, Y.; Ohashi, H.; Pettersson, L. G. M.; Nilsson, A.; Shin, S. Chem. Phys. Lett. 2008, 460, 387–400. (4) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Science 2004, 306, 851–853. (5) Head-Gordon, T.; Johnson, M. E. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 7973–7977. (6) Mantz, Y. A.; Chen, B.; Martyna, G. J. Chem. Phys. Lett. 2005, 405, 294–299.

IR Spectroscopy and H-Bond Dynamics of Liquid Water (7) Mantz, Y. A.; Chen, B.; Martyna, G. J. J. Phys. Chem. B 2006, 110, 3540–3554. (8) Soper, A. K. J. Phys.: Condens. Matter 2005, 17, S3273–S3282. (9) Leetmaa, M.; Ljungberg, M.; Ogasawara, H.; Odelius, M.; Na¨slund, L.-A.; Nilsson, A.; Pettersson, L. G. M. J. Chem. Phys. 2006, 125, 244510. (10) Prendergast, D.; Galli, G. Phys. ReV. Lett. 2006, 96, 215502. (11) Leetmaa, M.; Wikfeldt, K. T.; Ljungberg, M. P.; Odelius, M.; Swenson, J.; Nilsson, A.; Pettersson, L. G. M. J. Chem. Phys. 2008, 129, 084502. (12) Smith, J. D.; Cappa, C. D.; Wilson, K. R.; Messer, B. M.; Cohen, R. C.; Saykally, R. J. Science 2004, 306, 851–853. (13) Wang, Z. H.; Pakoulev, A.; Pang, Y.; Dlott, D. D. Chem. Phys. Lett. 2003, 378, 281–288. (14) 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–14174. (15) Geissler, P. L. J. Am. Chem. Soc. 2005, 127, 14930–14935. (16) Walrafen, G. E.; Fisher, M. R.; Hokmabadi, M. S.; Yang, W. H. J. Chem. Phys. 1986, 85, 6970–6982. (17) Hart, R. T.; Benmore, C. J.; Neuefeind, J. C.; Kohara, S.; Tomberli, B.; Egelstaff, P. A. Phys. ReV. Lett. 2005, 94, 047801. (18) 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. (19) Soper, A. K.; Benmore, C. J. Phys. ReV. Lett. 2008, 101, 065502. (20) Tauber, M. J.; Mathies, R. A. J. Am. Chem. Soc. 2003, 125, 1394– 1402. (21) Walrafen, G. E. In Water a comprehensiVe treatise; Franks, F., Ed.; Plenum Press: New York, 1972; Vol. 1, pp 151-214. (22) Chumaevskii, N. A.; Rodnikova, M. N.; Sirotkin, D. A. J. Mol. Liq. 1999, 82, 39–46. (23) Schmidt, D. A.; Miki, K. J. Phys. Chem. A 2007, 111, 10119– 10122. (24) Rey, R.; Moller, K. B.; Hynes, J. T. J. Phys. Chem. A 2002, 106, 11993–11996. (25) Auer, B. M.; Skinner, J. L. J. Chem. Phys. 2008, 128, 224511. (26) Tan, H. S.; Piletic, I. R.; Fayer, M. D. J. Chem. Phys. 2005, 122, 174501. (27) Asbury, J. B.; Steinel, T.; Fayer, M. D. J. Phys. Chem. B 2004, 108, 6544–6554. (28) Asbury, J. B.; Steinel, T.; Kwak, K.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. J. Chem. Phys. 2004, 121, 12431–12446. (29) 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– 1119. (30) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Fayer, M. D. J. Chem. Phys. 2003, 119, 12981–12997. (31) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Goun, A.; Fayer, M. D. Phys. ReV. Lett. 2003, 91, 237402. (32) 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–5300. (33) Moilanen, D. E.; Wong, D.; Rosenfeld, D. E.; Fenn, E. E.; Fayer, M. D. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 375–380. (34) Park, S.; Moilanen, D. E.; Fayer, M. D. J. Phys. Chem. B 2008, 112, 5279–5290. (35) Piletic, I. R.; Moilanen, D. E.; Levinger, N. E.; Fayer, M. D. J. Am. Chem. Soc. 2006, 128, 10366–10367. (36) Piletic, I. R.; Tan, H. S.; Fayer, M. D. J. Phys. Chem. B 2005, 109, 21273–21284. (37) Steinel, T.; Asbury, J. B.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. Chem. Phys. Lett. 2004, 386, 295–300. (38) Steinel, T.; Asbury, J. B.; Zheng, J. R.; Fayer, M. D. J. Phys. Chem. A 2004, 108, 10957–10964. (39) Tan, H. S.; Piletic, I. R.; Riter, R. E.; Levinger, N. E.; Fayer, M. D. Phys. ReV. Lett. 2005, 94, 057405. (40) Zheng, J. R.; Fayer, M. D. J. Am. Chem. Soc. 2007, 129, 4328– 4335. (41) Bakker, H. J. Chem. ReV. 2008, 108, 1456–1473. (42) Bakker, H. J.; Kropman, M. F.; Omta, A. W. J. Phys.: Condens. Matter 2005, 17, S3215–S3224. (43) Bakker, H. J.; Rezus, Y. L. A.; Timmer, R. L. A. J. Phys. Chem. A 2008, 112, 11523–11534. (44) Bakker, H. J.; Woutersen, S.; Nienhuys, H. K. Chem. Phys. 2000, 258, 233–245. (45) Dokter, A. M.; Woutersen, S.; Bakker, H. J. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 15355–15358. (46) Dokter, A. M.; Woutersen, S.; Bakker, H. J. J. Chem. Phys. 2007, 126, 124507. (47) Gilijamse, J. J.; Lock, A. J.; Bakker, H. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 3202–3207. (48) Kropman, M. F.; Bakker, H. J. Science 2001, 291, 2118–2120. (49) Kropman, M. F.; Bakker, H. J. Chem. Phys. Lett. 2003, 370, 741– 746.

J. Phys. Chem. B, Vol. 113, No. 39, 2009 13129 (50) Kropman, M. F.; Bakker, H. J. J. Am. Chem. Soc. 2004, 126, 9135– 9141. (51) Kropman, M. F.; Nienhuys, H. K.; Bakker, H. J. Phys. ReV. Lett. 2002, 88, 077601. (52) Kropman, M. F.; Nienhuys, H. K.; Woutersen, S.; Bakker, H. J. J. Phys. Chem. A 2001, 105, 4622–4626. (53) Lock, A. J.; Bakker, H. J. J. Chem. Phys. 2002, 117, 1708–1713. (54) Nienhuys, H. K.; Lock, A. J.; van Santen, R. A.; Bakker, H. J. J. Chem. Phys. 2002, 117, 8021–8029. (55) Nienhuys, H. K.; van Santen, R. A.; Bakker, H. J. J. Chem. Phys. 2000, 112, 8487–8494. (56) Nienhuys, H. K.; Woutersen, S.; van Santen, R. A.; Bakker, H. J. J. Chem. Phys. 1999, 111, 1494–1500. (57) Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J. J. Chem. Phys. 2003, 119, 12457–12461. (58) Omta, A. W.; Kropman, M. F.; Woutersen, S.; Bakker, H. J. Science 2003, 301, 347–349. (59) Rezus, Y. L. A.; Bakker, H. J. J. Chem. Phys. 2005, 123, 114502. (60) Rezus, Y. L. A.; Bakker, H. J. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 18417–18420. (61) Rezus, Y. L. A.; Bakker, H. J. Phys. ReV. Lett. 2007, 99, 148301. (62) Rezus, Y. L. A.; Bakker, H. J. J. Phys. Chem. A 2008, 112, 2355– 2361. (63) Tielrooij, K. J.; Petersen, C.; Rezus, Y. L. A.; Bakker, H. J. Chem. Phys. Lett. 2009, 471, 71–74. (64) Woutersen, S.; Bakker, H. J. Phys. ReV. Lett. 1999, 83, 2077– 2080. (65) 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–13022. (66) Eaves, J. D.; Tokmakoff, A.; Geissler, P. L. J. Phys. Chem. A 2005, 109, 9424–9436. (67) Fecko, C. J.; Eaves, J. D.; Loparo, J. J.; Tokmakoff, A.; Geissler, P. L. Science 2003, 301, 1698–1702. (68) Fecko, C. J.; Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2005, 122, 054506. (69) Loparo, J. J.; Fecko, C. J.; Eaves, J. D.; Roberts, S. T.; Tokmakoff, A. Phys. ReV. B 2004, 70, 180201. (70) Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 194521. (71) Loparo, J. J.; Roberts, S. T.; Tokmakoff, A. J. Chem. Phys. 2006, 125, 194522. (72) Nicodemus, R. A.; Tokmakoff, A. Chem. Phys. Lett. 2007, 449, 130–134. (73) Morrone, J. A.; Car, R. Phys. ReV. Lett. 2008, 101, 017801-4. (74) Silvestrelli, P. L.; Parrinello, M. Phys. ReV. Lett. 1999, 82, 3308. (75) Laasonen, K.; Sprik, M.; Parrinello, M. J. Chem. Phys. 1993, 99, 9080. (76) Fois, E. S.; Sprik, M.; Parrinello, M. Chem. Phys. Lett. 1994, 223, 411. (77) Sprik, M.; Hutter, J.; Parrinello, M. J. Chem. Phys. 1996, 195, 1142. (78) VandeVondele, J.; Mohamed, F.; Krack, M.; Hutter, J.; Sprik, M.; Parrinello, M. J. Chem. Phys. 2005, 122, 014515. (79) Silvestrelli, P. L.; Bernasconi, M.; Parrinello, M. Chem. Phys. Lett. 1997, 277, 478. (80) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 3572– 3580. (81) Raiteri, P.; Laio, A.; Parrinello, M. Phys. ReV. Lett. 2004, 93, 087801. (82) Sharma, M.; Resta, R.; Car, R. Phys. ReV. Lett. 2005, 95, 187401. (83) Chen, B.; Ivanov, I.; Klein, M. L.; Parrinello, M. Phys. ReV. Lett. 2003, 91, 215503. (84) Lee, H.-S.; Tuckerman, M. E. J. Phys. Chem. A 2006, 110, 5549– 5560. (85) Lee, H.-S.; Tuckerman, M. E. J. Chem. Phys. 2006, 125, 154507– 14. (86) Lee, H.-S.; Tuckerman, M. E. J. Chem. Phys. 2007, 126, 164501– 16. (87) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. J. Chem. Phys. 2004, 120, 300–311. (88) Izvekov, S.; Voth, G. A. J. Chem. Phys. 2002, 116, 10372–10376. (89) Morrone, J. A.; Srinivasan, V.; Sebastiani, D.; Car, R. J. Chem. Phys. 2007, 126, 234504–9. (90) Kuo, I.-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; Krack, M.; Parrinello, M. J. Phys. Chem. B 2004, 108, 12990–12998. (91) Kuo, I. F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I. J. Chem. Theory Comput. 2006, 2, 1274–1281. (92) Yoo, S.; Zeng, X. C.; Xantheas, S. S. J. Chem. Phys. 2009, 130, 221102. (93) Benighaus, T.; DiStasio, R. A.; Lochan, R. C.; Chai, J. D.; HeadGordon, M. J. Phys. Chem. A 2008, 112, 2702–2712.

13130

J. Phys. Chem. B, Vol. 113, No. 39, 2009

(94) Gianturco, F. A.; Paesani, F. In Conceptual perspectiVes in Quantum Chemistry; Calais, J.-L., Kryachko, E., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1997. (95) Gianturco, F. A.; Paesani, F. J. Chem. Phys. 2000, 113, 3011– 3019. (96) Gianturco, F. A.; Paesani, F. J. Chem. Phys. 2001, 115, 249–256. (97) Gianturco, F. A.; Paesani, F. Mol. Phys. 2001, 99, 689–698. (98) Gianturco, F. A.; Paesani, F.; Laranjeira, M. F.; Vassilenko, V.; Cunha, M. A. J. Chem. Phys. 1999, 110, 7832–7845. (99) Grimme, S. J. Comput. Chem. 2004, 25, 1463–1473. (100) Wu, X.; Vargas, M. C.; Nayak, S.; Lotrich, V.; Scoles, G. J. Chem. Phys. 2001, 115, 8748–8757. (101) Chai, J. D.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2008, 10, 6615–6620. (102) Lin, I. C.; Seitsonen, A. P.; Coutinho-Neto, M. D.; Tavernelli, I.; Rothlisberger, U. J. Phys. Chem. B 2009, 113, 1127–1131. (103) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2003, 118, 264– 272. (104) Corcelli, S. A.; Skinner, J. L. J. Phys. Chem. A 2005, 109, 6154– 6165. (105) Lawrence, C. P.; Skinner, J. L. Chem. Phys. Lett. 2003, 369, 472– 477. (106) Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2002, 117, 8847– 8854. (107) Moller, K. B.; Rey, R.; Hynes, J. T. J. Phys. Chem. A 2004, 108, 1275–1289. (108) Harder, E.; Eaves, J. D.; Tokmakoff, A.; Berne, B. J. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 11611–11616. (109) Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L. J. Chem. Phys. 2004, 120, 8107–8117. (110) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2004, 121, 8887–8896. (111) Hayashi, T.; Jansen, T. L.; Zhuang, W.; Mukamel, S. J. Phys. Chem. A 2005, 109, 64–82. (112) Jansen, T. L.; Hayashi, T.; Zhuang, W.; Mukamel, S. J. Chem. Phys. 2005, 123, 114504. (113) Laage, D.; Hynes, J. T. Chem. Phys. Lett. 2006, 433, 80–85. (114) Schmidt, J. R.; Roberts, S. T.; Loparo, J. J.; Tokmakoff, A.; Fayer, M. D.; Skinner, J. L. Chem. Phys. 2007, 341, 143–157. (115) Wu, Y. J.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (116) Fernandez, R. G.; Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2006, 124, 144506. (117) Auer, B.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 14215–14220. (118) Paesani, F.; Voth, G. A. J. Phys. Chem. B 2009, 113, 5702–5719. (119) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5093–5105. (120) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 100, 5106–5117. (121) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6157–6167. (122) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6168–6183. (123) Cao, J.; Voth, G. A. J. Chem. Phys. 1994, 101, 6184–6192. (124) Jang, S.; Voth, G. A. J. Chem. Phys. 1999, 111, 2357–2370. (125) Jang, S.; Voth, G. A. J. Chem. Phys. 1999, 111, 2371–2384. (126) Fanourgakis, G. S.; Xantheas, S. S. J. Chem. Phys. 2008, 128, 074506. (127) Lobaugh, J.; Voth, G. A. J. Chem. Phys. 1997, 106, 2400–2410. (128) Hone, T. D.; Voth, G. A. J. Chem. Phys. 2004, 121, 6412–6422. (129) Paesani, F.; Zhang, W.; Case, D. A.; Cheatham, T. E., III; Voth, G. A. J. Chem. Phys. 2006, 125, 184507. (130) Paesani, F.; Iuchi, S.; Voth, G. A. J. Chem. Phys. 2007, 127, 074506. (131) Paesani, F.; Voth, G. A. J. Phys. Chem. C 2008, 112, 324–327. (132) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 7146–7146. (133) Herna´ndez de la Pen˜a, L.; Kusalik, P. G. J. Chem. Phys. 2004, 121, 5992–6002.

Paesani et al. (134) Herna´ndez de la Pen˜a, L.; Kusalik, P. G. J. Am. Chem. Soc. 2005, 127, 5246–5251. (135) Herna´ndez de la Pen˜a, L.; Razul, M. S. G.; Kusalik, P. G. J. Phys. Chem. A 2005, 109, 7236–7241. (136) Herna´ndez de la Pen˜a, L.; Kusalik, P. G. J. Chem. Phys. 2006, 125, 054512. (137) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361– 9381. (138) Jang, S.; Jang, S.; Voth, G. A. J. Phys. Chem. A 1999, 103, 9512– 9520. (139) Habershon, S.; Fanourgakis, G. S.; Manolopoulos, D. E. J. Chem. Phys. 2008, 129, 074501. (140) Voth, G. A. AdV. Chem. Phys. 1996, 93, 135–218. (141) Reichman, D. R.; Roy, P.-N.; Jang, S.; Voth, G. A. J. Chem. Phys. 2000, 113, 919–929. (142) Krishna, V.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18953– 18957. (143) Jang, S. J. Chem. Phys. 2006, 124, 064107-8. (144) Paesani, F.; Voth, G. A. J. Chem. Phys. 2008, 129, 194113–10. (145) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1479– 1492. (146) Xantheas, S. S.; Burnham, C. J.; Harrison, R. J. J. Chem. Phys. 2002, 116, 1493–1499. (147) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1500– 1510. (148) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115– 5124. (149) Fanourgakis, G. S.; Xantheas, S. S. J. Phys. Chem. A 2006, 110, 4100–4106. (150) Thole, B. T. Chem. Phys. 1981, 59, 341–350. (151) Xantheas, S. S.; Apra`, E. J. Chem. Phys. 2004, 120, 823–828. (152) Fanourgakis, G. S.; Apra, E.; Xantheas, S. S. J. Chem. Phys. 2004, 121, 2655–2663. (153) Bulusu, S.; Yoo, S.; Apra, E.; Xantheas, S.; Zeng, X. C. J. Phys. Chem. A 2006, 110, 11781–11784. (154) Hankins, D.; Moskowitz, J. W.; Stillinger, F. H. J. Chem. Phys. 1970, 53, 4544–4554. (155) Xantheas, S. S. Chem. Phys. 2000, 258, 225–231. (156) Xantheas, S. S. J. Chem. Phys. 1994, 100, 7523–7534. (157) Partridge, H.; Schwenke, D. W. J. Chem. Phys. 1997, 106, 4618– 4639. (158) Fanourgakis, G. S.; Xantheas, S. S. J. Chem. Phys. 2006, 124, 174504. (159) Fanourgakis, G. S.; Schenter, G. K.; Xantheas, S. S. J. Chem. Phys. 2006, 125, 141102. (160) Numerov, B. V. Mon. Not. R. Astron. Soc. 1924, 84, 592. (161) Schenter, G. K.; Garrett, B. C.; Voth, G. A. J. Chem. Phys. 2000, 113, 5171–5178. (162) Bertie, J. E.; Lan, Z. Appl. Spectrosc. 1996, 50, 1047–1057. (163) Venyaminov, S. Y.; Prendergast, F. G. Anal. Biochem. 1997, 248, 234–245. (164) Guillot, B. J. Chem. Phys. 1991, 95, 1543–1551. (165) Skinner, J. L. Mol. Phys. 2008, 106, 2245–2253. (166) Gruetzmacher, J. A.; Nome, R. A.; Moran, A. M.; Scherer, N. F. J. Chem. Phys. 2008, 129, 224502–10. (167) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. J. Chem. Phys. 2005, 123, 044513. (168) Hermansson, K. J. Chem. Phys. 1993, 99, 861–868. (169) Corcelli, S. A.; Lawrence, C. P.; Asbury, J. B.; Steinel, T.; Fayer, M. D.; Skinner, J. L. J. Chem. Phys. 2004, 121, 8897–8900. (170) Smith, J. D.; Saykally, R. J.; Geissler, P. L. J. Am. Chem. Soc. 2007, 129, 13847–13856. (171) Yeremenko, S.; Pshenichnikov, M. S.; Wiersma, D. A. Chem. Phys. Lett. 2003, 369, 107–113.

JP907648Y