Mixed quantum/classical approach to OH-stretch inelastic incoherent neutron scattering spectroscopy for ambient and supercooled liquid water and ice Ih L. Shi, and J. L. Skinner
Citation: The Journal of Chemical Physics 143, 014503 (2015); doi: 10.1063/1.4923387 View online: http://dx.doi.org/10.1063/1.4923387 View Table of Contents: http://aip.scitation.org/toc/jcp/143/1 Published by the American Institute of Physics
Articles you may be interested in Evidence for a liquid-liquid critical point in supercooled water within the E3B3 model and a possible interpretation of the kink in the homogeneous nucleation line The Journal of Chemical Physics 144, 214501 (2016); 10.1063/1.4952991 IR spectra of water droplets in no man’s land and the location of the liquid-liquid critical point The Journal of Chemical Physics 145, 124509 (2016); 10.1063/1.4963736 IR and SFG vibrational spectroscopy of the water bend in the bulk liquid and at the liquid-vapor interface, respectively The Journal of Chemical Physics 143, 014502 (2015); 10.1063/1.4923462 Infrared and Raman line shapes for ice Ih. II. H2O and D2O The Journal of Chemical Physics 133, 244504 (2010); 10.1063/1.3516460 IR and Raman spectra of liquid water: Theory and interpretation The Journal of Chemical Physics 128, 224511 (2008); 10.1063/1.2925258 Infrared and Raman line shapes for ice Ih. I. Dilute HOD in and The Journal of Chemical Physics 132, 204505 (2010); 10.1063/1.3430518
THE JOURNAL OF CHEMICAL PHYSICS 143, 014503 (2015)
Mixed quantum/classical approach to OH-stretch inelastic incoherent neutron scattering spectroscopy for ambient and supercooled liquid water and ice Ih L. Shi and J. L. Skinner Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, USA
(Received 27 April 2015; accepted 22 June 2015; published online 7 July 2015) OH-stretch inelastic incoherent neutron scattering (IINS) has been measured to determine the vibrational density of states (VDOS) in the OH-stretch region for liquid water, supercooled water, and ice Ih, providing complementary information to IR and Raman spectroscopies about hydrogen bonding in these phases. In this work, we extend the combined electronic-structure/molecular-dynamics (ES/MD) method, originally developed by Skinner and co-workers to simulate OH-stretch IR and Raman spectra, to the calculation of IINS spectra with small k values. The agreement between theory and experiment in the limit k → 0 is reasonable, further validating the reliability of the ES/MD method in simulating OH-stretch spectroscopy in condensed phases. The connections and differences between IINS and IR spectra are analyzed to illustrate the advantages of IINS over IR in estimating the OH-stretch VDOS. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4923387] I. INTRODUCTION
Understanding water in its condensed phases (liquid water, crystalline, and amorphous ices) has been a recurring research theme for over a century. Vibrational spectroscopy, both of the low-frequency intermolecular modes and the highfrequency intramolecular modes, has played an important role in this understanding. In this paper, we discuss vibrational spectroscopy in the OH-stretch region, which provides useful information about hydrogen bonding because the OH-stretch frequency is very sensitive to local environment. This field has been dominated by IR and Raman spectroscopy. For example, the IR spectrum of liquid water shows a peak at about 3400 cm−1 (and also has a small shoulder at about 3250 cm−1) with a full-width-half-maximum (FWHM) of about 375 cm−1,1 manifesting the effect of hydrogen bonding on the OH-stretch frequency, and the structural heterogeneity in liquid water.2 The interpretation of this spectrum is not straightforward, due in part to the presence of strong intermolecular vibrational coupling.2,3 Inelastic incoherent neutron scattering (IINS) is an alternative vibrational spectroscopy. This technique is complementary to IR and Raman, and provides more direct information about the vibrational density of states (VDOS). Due to the availability of high-energy neutron sources, IINS now can span the entire infrared region for typical molecular (fundamental) vibrations.4–6 Since the 1980s, OH-stretch IINS has been measured for liquid water, ice Ih, supercooled water, amorphous ice, and supercritical water in hopes of shedding some light on the hydrogen bonding in these phases, as well as on the interpretations of the complex IR and Raman spectra.7–16 In fact, the OH-stretch IINS results are different from their optical counterparts for these phases, and a clear theoretical explanation for these differences still seems lacking. 0021-9606/2015/143(1)/014503/6/$30.00
Considering the rich information in the OH-stretch IINS of water, a direct simulation of such is desired. However, just as in the case of OH-stretch IR and Raman spectra, IINS requires a quantum-mechanical treatment of the OH stretches as their vibrational frequencies (over 3000 cm−1) are much larger than thermal energy at most relevant temperatures (at 300 K, k BT is about 200 cm−1). There are approaches in the literature to calculate IINS semiclassically,17–20 but, to the best of our knowledge, they have not been applied to the OH-stretch IINS of water. In the 1990s, Bratos and coworkers have put forth a mixed quantum/classical method to calculate the OH-stretch IINS spectra for water,21–24 in which some empirical parameters were needed in the calculations and intermolecular OH-stretch vibrational coupling was not considered explicitly. Quite recently Bowman and co-workers have calculated, fully quantum mechanically, the VDOS up to 4000 cm−1 using their very accurate WHBB potential surface.25 In recent years, Skinner and co-workers have developed a combined electronic-structure/molecular-dynamics (ES/MD) method to simulate OH-stretch vibrational spectroscopy,26,27 and it has been successful in reproducing various OH (OD) stretch experimental spectra (e.g., IR, Raman, 2DIR, pumpprobe anisotropy, sum-frequency generation) of water in many different molecular environments (e.g., liquid water, ice Ih, air/water interface, hexamer and amorphous ices).2,26,28–43 The goal of the present work is to extend this ES/MD method to the simulation of the OH-stretch IINS of water. The rest of the paper is organized as follows. In Sec. II, we outline the general theory of IINS and the ES/MD method for IINS with small k values. In Sec. III, the calculated IINS spectra in the limit k → 0 for liquid water, supercooled water, and ice Ih are compared with experiment. In Sec. IV, the connections and differences between IR and IINS spectra are elucidated in detail. In Sec. V, we conclude.
143, 014503-1
© 2015 AIP Publishing LLC
014503-2
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
II. MIXED QUANTUM/CLASSICAL EXPRESSION FOR THE IINS SPECTRUM
For H-containing systems like water, since the IINS crosssection of hydrogen is much larger than those of other atoms including deuterium,5 it is a rather good approximation to consider only the scattering due to hydrogen in the IINS calculation. Assuming the first Born approximation and a Fermi pseudopotential between the incident neutron and the scattering atom, the IINS double-differential cross-section is proportional to the incoherent scattering function (sometimes also called the incoherent dynamical structure factor), given by44 ∞ 1 dt e−iωt I(⃗k,t), (1) S(⃗k,ω) = 2π −∞ where ⃗k is the momentum transferred during the scattering, and I(⃗k,t) is the incoherent intermediate scattering function, defined by44 1 −i k⃗ ·⃗r l (0) i k⃗ ·⃗r l (t) e , (2) I(⃗k,t) = e N l where the brackets denote an equilibrium quantum-mechanical ensemble average, N is the number of hydrogen atoms, and ⃗r l (t) is the time-dependent position operator associated with the lth hydrogen in the Heisenberg picture. Experimentally one often measures the scattering for low k (≡ |⃗k|), and defines the quantity ω2 ⃗ S( k,ω), k→ 0 k 2
G0(ω) ≡ lim
(3)
which can be called the IINS spectrum. As stated in the Introduction, it is essential to treat the OH stretch quantum-mechanically, and the ES/MD method developed by Skinner and co-workers is an accurate and efficient mixed quantum/classical approach to do so.27 In this method, the quantum-mechanical system consists of many coupled local OH stretches (also called chromophores), and its many-body vibrational ground state is denoted as |0⟩. When only the jth OH chromophore is excited to its first vibrational excited state, the vibrational state of the system is | j⟩ (note that higher excited states may be reached in nonlinear spectroscopy). The Hamiltonian (up to and including only the oneexciton band) is therefore H = H0|0⟩⟨0| + Hi j |i⟩ ⟨ j|. (4) ij
Note that in our treatment we do not include bend degrees of freedom. The system also has low-frequency modes (the bath, in this case translations and rotations of water molecules) and H0 and Hi j are operators in the bath Hilbert space. Using the completeness relation (for the ground vibrational state and the one-exciton band), 1 = |0⟩ ⟨0| + i |i⟩⟨i|, this can be rewritten as H = H0 + ~ κ i j |i⟩ ⟨ j|, (5) ij
where ~κ i j = (Hii − H0)δ i j + Hi j (1 − δ i j ).
(6)
Thus, the diagonal elements are the (typically anharmonic) local-mode transition frequencies, and the off-diagonal elements are the vibrational couplings (both intra- and intermolecular). Both frequencies and couplings are formally bath operators. H0 is the bath Hamiltonian. In what follows, the bath will become classical, and so the frequencies and couplings become functions of the bath coordinates. For a given bath configuration, these frequencies and couplings are evaluated through so-called spectroscopic “maps” based on extensive electronic-structure calculations.45 The bath coordinates evolve through a classical MD simulation, and so the frequencies and couplings fluctuate in time. In the Appendix we show that within this mixed quantum/classical approach, the quantum time-correlation function of two operators A and B can be written approximately as ⟨A0 j (0)Fj i (t)Bi0(t)⟩, (7) ⟨A(0)B(t)⟩ = ij
where now on the r.h.s. the brackets indicate a classical average, the matrix elements A0 j and Bi0 are bath variables, and the matrix F(t) satisfies the equation ˙ = iF(t)κ(t), F(t)
(8)
subject to the initial condition that Fi j (0) = δ i j . (κ i j are also bath variables.) The intermediate scattering function is therefore given approximately by 1 ⃗ ⃗ ⟨0|e−i k ·⃗r l (0)| j⟩Fj i (t)⟨i|ei k ·⃗r l (t)|0⟩ . I(⃗k,t) = N i jl
(9)
For the instantaneous position operator of hydrogen atom l, we then write ⃗l + uˆl x l , ⃗r l = R
(10)
where x l is the displacement from the equilibrium value of the H atom along the OH bond, uˆl is the unit vector of the OH ⃗l (t) is the equilibrium position of the H atom. To bond, and R first order in k, ⃗
l ⟨0|e−i k ·⃗r l | j⟩ = −ikδ jl kˆ · uˆl x 01 ,
(11)
l where x 01 is the 0-1 matrix element of the displacement coordinate for hydrogen l. For an isotropic system, G0(ω) then becomes ∞ ω2 dt e−iωt G0(ω) = 6N π −∞
l l × uˆl (0) · uˆl (t)x 01 (0)x 01 (t)Fll (t) e−|t |/2T1. (12) l
We have included the effect of the vibrational lifetime, T1, phenomenologically.27 Thus, in the end, the IINS spectrum is related to a classical time-correlation function involving the fluctuating frequencies and vibrational couplings, the fluctuating direction of each OH bond, and the fluctuating position matrix elements (which depend only weakly on the local environment).
014503-3
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
III. IINS SPECTRA FOR WATER AND ICE
To perform theoretical IINS calculations, the time evolution of the bath degrees of freedom was carried out by classical MD simulation with the third version of the E3B water model (E3B3),46 which uses the TIP4P/2005 model47 as a reference potential, and includes three-body interactions explicitly.48,49 E3B3 shows improved performance in describing some properties of water over the previous version (E3B2): for instance, better melting temperature of ice (260 K for E3B3 versus 251 K for E3B2, compared to the experimental value of 273 K), and better ice density at one atmosphere and 273 K (0.918 g/ml for E3B3 versus 0.935 g/ml for E3B2, compared to the experimental value of 0.917 g/ml).46 All simulations of liquid water, supercooled water, and ice Ih were performed in the NPT ensemble with 432-molecule simulation boxes using a modified GROMACS package 4.5.5.50–53 The pressure was maintained at 1 bar using a Parrinello-Rahman barostat54 with a coupling constant of 10.0 ps in production runs, and the temperatures were controlled by a Nosé-Hoover thermostat55,56 with a coupling constant of 5.0 ps. The geometry of the water molecule was constrained using the SETTLE algorithm,57 and the electrostatic interactions with periodic boundary conditions were treated using the particle-mesh Ewald (PME) method.58,59 The cutoff for the Lennard-Jones potential between oxygen atoms was 1.0 nm, and the longrange dispersion corrections were applied to the energy and pressure.60 A proton-disordered configuration generated by Hayward and Reimers61 was employed as the initial configuration of ice Ih. The simulation time step was 1 fs, and the trajectories were saved every 5 fs for spectral calculations. The spectroscopic maps used in the present study are from Table 1 of Ref. 45, and in order to use these TIP4P-based maps the geometry of the E3B3 (TIP4P/2005) water was adjusted to be that of TIP4P62 for the purpose of the spectral calculations. Fig. 1 displays the calculated G0(ω) (from Eq. (12)) for liquid water (top panel), supercooled water (middle panel), and (poly-crystalline) ice Ih (bottom panel), along with extrapolated results from experiment for the same systems. In the calculations, for the vibrational lifetimes we used 260 fs for ambient and supercooled water,63 and 300 fs for ice Ih.64 Note that as the melting temperature of the E3B3 model is 260 K,46 we choose 258 K for ice Ih and supercooled water in our simulations. All the theoretical and experimental spectra show essentially single peaks, except that the theoretical spectrum for ice Ih shows a shoulder on the low-frequency side of the main peak. This shoulder may also exist in experiment (it does actually exist in the 15 K experiments by Li10 but these results are not extrapolated to k = 0), but the modest frequency resolution of the experiment precludes a definitive answer. In terms of peak position, fair agreement between theory and experiment is seen for liquid and supercooled water. However, for ice Ih, the calculated spectrum is blue-shifted by about 50 cm−1 compared to the experiment. This blue shift is probably due to the fact that the frequency map overestimates the OH-stretch frequency in ice Ih by about 30-40 cm−1.45 Regarding the linewidth, our theory is somewhat too narrow for supercooled water and ice Ih, and this discrepancy might be due to the inadequacy of the experimental extrapolation, the experimental resolution, or the
FIG. 1. Experimentally extrapolated (black lines) and calculated (red lines) G0(ω) for OH stretch in liquid water (top panel), supercooled water (middle panel), and ice Ih (bottom panel). Reference numbers for experiments are: (a) Ref. 16, (b) Ref. 11, (c) Ref. 14. All the spectra are normalized to have the same height of 1.
fact that our simulation temperatures (258 K) are lower than corresponding experimental ones (∼270 K). Overall, however, the agreement between theory and experiment is reasonable, especially in reproducing the red peak-shifts from liquid water, to supercooled water, to ice Ih (from top to bottom). IV. COMPARISON OF IINS AND IR SPECTRA
OH-stretch IR and Raman spectroscopies have been used to study liquid water and ice for years, and there are a number of experimental spectra available in the literature (see citations in Refs. 2 and 34). From the top and bottom panels of Fig. 1, one might already realize one difference between experimental IR line-shapes and G0(ω) for water: the IR spectra are peaked at lower frequencies (about 3400 cm−1 for liquid water at room temperature,1 and about 3250 cm−1, for ice at about 270 K34,65) compared to the corresponding G0(ω) (about 3450 cm−1 for liquid water,12,16 and about 3350 cm−1 for ice14,16). In this section, we will review the theoretical results for IR spectra, within the same mixed quantum/classical approach, and then highlight the two major differences between IR and IINS spectra. The IR spectrum for an isotropic system is related to the Fourier transform of the dipole-dipole time-correlation function,44 ∞ 1 ⃗ ⃗ I(ω) = dt e−iωt M(0) · M(t) , (13) 2π −∞ ⃗ is the total dipole of the sample (in the ground elecwhere M tronic state). The total dipole can be expanded to first order in the local-mode displacements. Within the bond-dipole approximation the required matrix element (according to Eq. (7)) is therefore ⃗ | j⟩ = µ′j uˆ j x j , ⟨0| M 01
(14)
014503-4
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
where µ′ is the dipole derivative, and the other symbols were defined previously. The spectrum, then, becomes2,3 ∞ 1 j i I(ω) = µ′j (0)µ′i (t)uˆ j (0) · uˆi (t)x 01 (0)x 01 (t)Fj i (t) e−|t |/2T1. dt e−iωt 2π −∞ ij
(15)
This expression involves a double sum over chromophores, and so there are two distinct contributions, when i = j and when i , j. Explicitly, they are ∞ 1 j j −iωt ′ ′ Iinc(ω) = dt e µ j (0)µ j (t)uˆ j (0) · uˆ j (t)x 01(0)x 01(t)Fj j (t) e−|t |/2T1 (16) 2π −∞ j and 1 Icoh(ω) = 2π
∞
dt e −∞
−iωt
µ′j (0)µ′i (t)uˆ j (0)
·
j i uˆi (t)x 01 (0)x 01 (t)Fj i (t)
e−|t |/2T1.
(17)
i, j
The first term is closely related to Eq. (12) for IINS, and so we call it “incoherent.” We call the second term “coherent” for two reasons: (1) it is similar to the major contribution to coherent neutron scattering; (2) this decomposition has been used previously to analyze IR line shapes,3 and this second term is responsible for “coherent” vibrational energy transfer. In Fig. 2, the calculated IR spectra (solid black lines) and their incoherent and coherent contributions (dashed and dotted black lines, respectively) for supercooled water (right panel) and ice Ih (left panel) at 258 K are shown in the top panels. One sees that the full line shape arises from a significant cancellation of the two terms. The coherent term would be zero in the absence of vibrational coupling and energy transfer.3 In the same figure, bottom panels, we show the calculated corresponding G0(ω) (red solid lines). It is evident that in each case the incoherent IR line shape Iinc(ω) closely resembles G0(ω). In fact, Eq. (12) for IINS and Eq. (16) for the incoherent contribution to IR are the same except for the factors of the dipole derivative, µ′, in the latter. µ′ is known to depend sensitively on the molecular environment, and stronger hydrogen bonding leads to a larger value.66 This “non-Condon” effect
has been shown to be important for various kinds of IR spectroscopy in liquid water. Since stronger hydrogen bonding corresponds to lower frequency, non-Condon effects tend to shift the spectrum to the red. Indeed, in Fig. 2 the incoherent contribution to the IR spectrum is red-shifted (slightly, by about 20 cm−1) compared to the IINS spectrum for supercooled water. In ice, however, all molecules have roughly the same hydrogen-bonding environment (they all have four hydrogen bonds), and so the dipole derivative is more or less constant. In this case, G0(ω) ∼ ω2 Iinc(ω), and indeed, in Fig. 2 we see that is the case (the very weak ω2 dependence can be ignored). To illustrate the relation between the VDOS and G0(ω), we consider the system in the inhomogeneous limit, when bath dynamics becomes unimportant.2,67 We also take T1 → ∞. In this limit, each of the fluctuating quantities at time t can be replaced by its value at time 0, and so, for example, we have ∞ l G0(ω)/ω2 ∼ dt e−iωt (0)2[eiκ(0)t ]ll . (18) x 01 −∞
l
FIG. 2. Calculated OH stretch IR and IINS spectra for ice Ih (left panel) and supercooled water (right panel) at 258 K. IR line shapes (solid black lines) with incoherent (dashed black lines, from Eq. (16)) and coherent (dotted black lines, from Eq. (17)) contributions are shown in the top, and IINS spectra G0(ω) (red lines, from Eq. (12)) are in the bottom. Frequency distributions (solid blue lines), calculated from Eq. (21), are also shown in the bottom panel.
014503-5
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
The matrix elements fluctuate only weakly, and therefore to an excellent approximation can be taken outside the average, in which case, ∞ 2 G0(ω)/ω ∼ dt e−iωt Tr[eiκ(0)t ] . (19) −∞
Since this result is now independent of the time origin, the argument of κ can be dropped. Further, κ can be diagonalized by an orthogonal matrix, QT κQ = λ,
(20)
where λ is diagonal, and the eigenvalues λ l are the frequencies of the vibrational eigenmodes. Therefore, it is easy to see that in this limit the IINS spectrum is 2 G0(ω)/ω ∼ δ(ω − λ l ) ≡ f (ω). (21) l
f (ω) is the VDOS. Note that this is not simply a harmonic result; rather, this is the VDOS even for arbitrarily large localmode anharmonicity. The calculated f (ω) for supercooled water and ice Ih at 258 K are shown as solid blue lines in Fig. 2, and the similarity between G0(ω) and f (ω) demonstrates again that IINS is useful in extracting the VDOS, although G0(ω) is sharpened from f (ω) due to motional narrowing. In contrast, the total IR line shape (solid black lines) is quite different from the VDOS both in terms of peak position and linewidth due to non-Condon effects and the coherent contribution.
ACKNOWLEDGMENTS
This work was supported by NSF Grant No. CHE1058752.
APPENDIX: MIXED QUANTUM/CLASSICAL EXPRESSION FOR QUANTUM TIME-CORRELATION FUNCTIONS
The equilibrium quantum time-correlation function for two arbitrary operators A and B can be written as ⟨A(0)B(t)⟩ = Tr[e−β H Aei H t /~ Be−i H t /~]/Tr[e−β H ],
(A1)
where in this case the Hamiltonian is given by Eq. (5) (and β = 1/k BT). The trace is over the bath states, and over the ground and one-exciton states. For these high-frequency vibrations ~κ ii ≫ k BT, and so the vibrational part of the trace involves only the ground state. Moreover, for the frequencies (in the eventual Fourier transform) of interest ω ∼ κ ii , this becomes A0 j ⟨ j|ei H t /~|i⟩Bi0e−i H0t /~]/ ⟨A(0)B(t)⟩ = Tr b [e−β H0 ij
Tr b [e
−β H0
],
(A2)
where Tr b indicates a trace over bath states. Bi0 is a bath operator, and so its time dependence is Bi0(t) = ei H0t /~ Bi0e−i H0t /~. With this definition, this becomes ⟨A(0)B(t)⟩ = Tr b [e−β H0 A0 j Fj i (t)Bi0(t)]/Tr b [e−β H0], ij
(A3) where
V. CONCLUDING REMARKS
In recent years, due to the availability of high-energy neutron sources, OH-stretch IINS spectra for water in condensed phases have been measured, and OH-stretch VDOS have been extracted from these experimental measurements, providing a good basis to evaluate modeling methods of OHstretch vibrational spectroscopy. In this work, we extended the ES/MD method, developed by Skinner and co-workers for optical vibrational spectroscopy (e.g., IR and Raman), to simulating OH-stretch IINS spectra at small k values. The reasonable agreement between theory and experiment for liquid water, supercooled water, and ice Ih further validates the ES/MD method, which has been quite successful in modeling various optical vibrational spectroscopies. We then go on to analyze the differences between IINS and IR spectra. We find there are two: first, the IR spectrum involves the dipole derivative, which tends to weight lower frequencies more; second and more importantly, the IR spectrum involves two terms—one that is analogous to incoherent neutron scattering and another that is analogous to coherent scattering. This second term involves the correlations and couplings between different molecules, and is signed (can be positive or negative). Adding this term to the “incoherent” contribution can dramatically change the spectrum. For both of these reasons the IR spectrum is more difficult to interpret than the IINS spectrum, which to a reasonably good approximation is just the VDOS.
Fj i (t) = ⟨ j|ei H t /~e−i H0t /~|i⟩.
(A4)
From Eq. (5), we have F(t) = ei(H0+~κ)t /~e−i H0t /~.
(A5)
From this it is clear that ˙ = iF(t)κ(t), F(t)
(A6)
where κ(t) = ei H0t /~ κe−i H0t /~. At this point we take the classical limit, such that F(t) and κ(t) are matrices in the oneexciton Hilbert space, but the matrix elements themselves are classical variables. In this limit the trace over the bath, weighted by the bath density operator, can be replaced by the classical equilibrium statistical mechanical average, recovering Eq. (7). 1J.
E. Bertie and Z. Lan, Appl. Spectrosc. 50, 1047 (1996). M. Auer and J. L. Skinner, J. Chem. Phys. 128, 224511 (2008). 3M. Yang and J. L. Skinner, Phys. Chem. Chem. Phys. 12, 982 (2010). 4P. C. H. Mitchell, S. F. Parker, A. J. Ramirez-Cuesta, and J. Tomkinson, Vibrational Spectroscopy with Neutrons With Applications in Chemistry, Biology, Materials Science and Catalysis (World Scientific Publishing Co. Pte. Ltd., Singapore, 2005). 5S. F. Parker, D. Lennon, and P. W. Albers, Appl. Spectrosc. 65, 1325 (2011). 6C. M. Brown, A. J. Ramirez-Cuesta, M. R. Johnson, and V. García-Sakai, Chem. Phys. 427, 1 (2013). 7S.-H. Chen, K. Toukan, C.-K. Loong, and D. L. Price, Phys. Rev. Lett. 53, 1360 (1984). 8C. Andreani, P. Bosi, F. Sacchetti, and C. K. Loong, J. Chem. Phys. 83, 750 (1985). 2B.
014503-6 9K.
L. Shi and J. L. Skinner
Toukan, M. A. Ricci, S.-H. Chen, C.-K. Loong, D. L. Price, and J. Teixeira, Phys. Rev. A 37, 2580 (1988). 10J.-C. Li, J. Chem. Phys. 105, 6733 (1996). 11M. A. Ricci, M. Nardone, A. Fontana, C. Andreani, and W. Hahn, J. Chem. Phys. 108, 450 (1998). 12D. D. Klug, C. A. Tulk, E. C. Svensson, and C.-K. Loong, Phys. Rev. Lett. 83, 2584 (1999). 13Y. Wang, A. I. Kolesnikov, S. L. Dong, and J. C. Li, Can. J. Phys. 81, 401 (2003). 14E. Perelli-Cippo, G. Gorini, M. Tardocchi, C. Andreani, A. Pietropaolo, R. Senesi, N. J. Rhodes, and E. M. Schooneveld, Appl. Phys. A 83, 453 (2006). 15C. Andreani, G. Romanelli, and R. Senesi, Chem. Phys. 427, 106 (2013). 16R. Senesi, D. Flammini, A. I. Kolesnikov, E. D. Murray, G. Galli, and C. Andreani, J. Chem. Phys. 139, 074504 (2013). 17E. J. Heller, J. R. Reimers, and G. Drolshagen, Phys. Rev. A 36, 2613 (1987). 18T. D. Hone and G. A. Voth, J. Chem. Phys. 121, 6412 (2004). 19I. R. Craig and D. E. Manolopoulos, Chem. Phys. 322, 236 (2006). 20J. Liu and W. H. Miller, J. Chem. Phys. 128, 144511 (2008). 21S. Bratos, G. Tarjus, M. Diraison, and J.-C. Leicknam, Phys. Rev. A 44, 2745 (1991). 22S. Bratos, M. Diraison, G. Tarjus, and J.-C. Leicknam, Phys. Rev. A 45, 5556 (1992). 23G. Tarjus, M. Diraison, J.-C. Leicknam, and S. Bratos, J. Mol. Struct. 296, 199 (1993). 24M. Diraison, J.-C. Leicknam, G. Tarjus, and S. Bratos, Phys. Rev. E 50, 2689 (1994). 25H. Liu, Y. Wang, and J. M. Bowman, J. Phys. Chem. B 117, 10046 (2013). 26S. A. Corcelli, C. P. Lawrence, and J. L. Skinner, J. Chem. Phys. 120, 8107 (2004). 27J. L. Skinner, B. M. Auer, and Y.-S. Lin, Adv. Chem. Phys. 142, 59 (2009). 28S. A. Corcelli and J. L. Skinner, J. Phys. Chem. A 109, 6154 (2005). 29B. M. Auer, R. Kumar, J. R. Schmidt, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 104, 14215 (2007). 30Y. Ni and J. L. Skinner, J. Chem. Phys. 141, 024509 (2014). 31F. Li and J. L. Skinner, J. Chem. Phys. 132, 204505 (2010). 32F. Li and J. L. Skinner, J. Chem. Phys. 133, 244504 (2010). 33F. Li and J. L. Skinner, J. Chem. Phys. 134, 099901 (2011). 34L. Shi, S. M. Gruenbaum, and J. L. Skinner, J. Phys. Chem. B 116, 13821 (2012). 35L. Shi and J. L. Skinner, J. Phys. Chem. B 117, 15536 (2013). 36L. Shi, F. Li, and J. L. Skinner, J. Chem. Phys. 140, 244503 (2014). 37I. V. Stiopkin, C. Weeraman, P. A. Pieniazek, F. Y. Shalhout, J. L. Skinner, and A. V. Benderskii, Nature 474, 192 (2011). 38P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Am. Chem. Soc. 133, 10360 (2011).
J. Chem. Phys. 143, 014503 (2015) 39P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Chem. Phys. 135, 044701
(2011). Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 110, 1992 (2013). 41C. J. Tainter, Y. Ni, L. Shi, and J. L. Skinner, J. Phys. Chem. Lett. 4, 12 (2013). 42C. J. Tainter, L. Shi, and J. L. Skinner, J. Chem. Phys. 140, 134503 (2014). 43T. L. C. Jansen, B. M. Auer, M. Yang, and J. L. Skinner, J. Chem. Phys. 132, 224503 (2010). 44D. A. McQuarrie, Statistical Mechanics (Harper and Row, New York, 1976). 45S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, J. Chem. Theory Comput. 9, 3109 (2013). 46C. J. Tainter, L. Shi, and J. L. Skinner, J. Chem. Theory Comput. 11, 2268 (2015). 47J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005). 48R. Kumar and J. L. Skinner, J. Phys. Chem. B 112, 8311 (2008). 49C. J. Tainter, P. A. Pieniazek, Y.-S. Lin, and J. L. Skinner, J. Chem. Phys. 134, 184501 (2011). 50B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). 51D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen, J. Comput. Chem. 26, 1701 (2005). 52E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Model. 7, 306 (2001). 53H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comput. Phys. Commun. 91, 43 (1995). 54M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982). 55S. Nosé, Mol. Phys. 52, 255 (1984). 56W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 57S. Miyamoto and P. A. Kollman, J. Comput. Chem. 13, 952 (1992). 58T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 59U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). 60M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1987). 61J. Hayward and J. Reimers, J. Chem. Phys. 106, 1518 (1997). 62W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 63A. J. Lock and H. J. Bakker, J. Chem. Phys. 117, 1708 (2002). 64R. L. A. Timmer and H. J. Bakker, J. Phys. Chem. A 114, 4148 (2010). 65S. G. Warren and R. E. Brandt, J. Geophys. Res. 113, D14220 doi:10.1029/2007jd009744 (2008). 66J. R. Schmidt, S. A. Corcelli, and J. L. Skinner, J. Chem. Phys. 123, 044513 (2005). 67B. M. Auer and J. L. Skinner, J. Chem. Phys. 127, 104105 (2007). 40Y.
Citation: The Journal of Chemical Physics 143, 014503 (2015); doi: 10.1063/1.4923387 View online: http://dx.doi.org/10.1063/1.4923387 View Table of Contents: http://aip.scitation.org/toc/jcp/143/1 Published by the American Institute of Physics
Articles you may be interested in Evidence for a liquid-liquid critical point in supercooled water within the E3B3 model and a possible interpretation of the kink in the homogeneous nucleation line The Journal of Chemical Physics 144, 214501 (2016); 10.1063/1.4952991 IR spectra of water droplets in no man’s land and the location of the liquid-liquid critical point The Journal of Chemical Physics 145, 124509 (2016); 10.1063/1.4963736 IR and SFG vibrational spectroscopy of the water bend in the bulk liquid and at the liquid-vapor interface, respectively The Journal of Chemical Physics 143, 014502 (2015); 10.1063/1.4923462 Infrared and Raman line shapes for ice Ih. II. H2O and D2O The Journal of Chemical Physics 133, 244504 (2010); 10.1063/1.3516460 IR and Raman spectra of liquid water: Theory and interpretation The Journal of Chemical Physics 128, 224511 (2008); 10.1063/1.2925258 Infrared and Raman line shapes for ice Ih. I. Dilute HOD in and The Journal of Chemical Physics 132, 204505 (2010); 10.1063/1.3430518
THE JOURNAL OF CHEMICAL PHYSICS 143, 014503 (2015)
Mixed quantum/classical approach to OH-stretch inelastic incoherent neutron scattering spectroscopy for ambient and supercooled liquid water and ice Ih L. Shi and J. L. Skinner Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, USA
(Received 27 April 2015; accepted 22 June 2015; published online 7 July 2015) OH-stretch inelastic incoherent neutron scattering (IINS) has been measured to determine the vibrational density of states (VDOS) in the OH-stretch region for liquid water, supercooled water, and ice Ih, providing complementary information to IR and Raman spectroscopies about hydrogen bonding in these phases. In this work, we extend the combined electronic-structure/molecular-dynamics (ES/MD) method, originally developed by Skinner and co-workers to simulate OH-stretch IR and Raman spectra, to the calculation of IINS spectra with small k values. The agreement between theory and experiment in the limit k → 0 is reasonable, further validating the reliability of the ES/MD method in simulating OH-stretch spectroscopy in condensed phases. The connections and differences between IINS and IR spectra are analyzed to illustrate the advantages of IINS over IR in estimating the OH-stretch VDOS. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4923387] I. INTRODUCTION
Understanding water in its condensed phases (liquid water, crystalline, and amorphous ices) has been a recurring research theme for over a century. Vibrational spectroscopy, both of the low-frequency intermolecular modes and the highfrequency intramolecular modes, has played an important role in this understanding. In this paper, we discuss vibrational spectroscopy in the OH-stretch region, which provides useful information about hydrogen bonding because the OH-stretch frequency is very sensitive to local environment. This field has been dominated by IR and Raman spectroscopy. For example, the IR spectrum of liquid water shows a peak at about 3400 cm−1 (and also has a small shoulder at about 3250 cm−1) with a full-width-half-maximum (FWHM) of about 375 cm−1,1 manifesting the effect of hydrogen bonding on the OH-stretch frequency, and the structural heterogeneity in liquid water.2 The interpretation of this spectrum is not straightforward, due in part to the presence of strong intermolecular vibrational coupling.2,3 Inelastic incoherent neutron scattering (IINS) is an alternative vibrational spectroscopy. This technique is complementary to IR and Raman, and provides more direct information about the vibrational density of states (VDOS). Due to the availability of high-energy neutron sources, IINS now can span the entire infrared region for typical molecular (fundamental) vibrations.4–6 Since the 1980s, OH-stretch IINS has been measured for liquid water, ice Ih, supercooled water, amorphous ice, and supercritical water in hopes of shedding some light on the hydrogen bonding in these phases, as well as on the interpretations of the complex IR and Raman spectra.7–16 In fact, the OH-stretch IINS results are different from their optical counterparts for these phases, and a clear theoretical explanation for these differences still seems lacking. 0021-9606/2015/143(1)/014503/6/$30.00
Considering the rich information in the OH-stretch IINS of water, a direct simulation of such is desired. However, just as in the case of OH-stretch IR and Raman spectra, IINS requires a quantum-mechanical treatment of the OH stretches as their vibrational frequencies (over 3000 cm−1) are much larger than thermal energy at most relevant temperatures (at 300 K, k BT is about 200 cm−1). There are approaches in the literature to calculate IINS semiclassically,17–20 but, to the best of our knowledge, they have not been applied to the OH-stretch IINS of water. In the 1990s, Bratos and coworkers have put forth a mixed quantum/classical method to calculate the OH-stretch IINS spectra for water,21–24 in which some empirical parameters were needed in the calculations and intermolecular OH-stretch vibrational coupling was not considered explicitly. Quite recently Bowman and co-workers have calculated, fully quantum mechanically, the VDOS up to 4000 cm−1 using their very accurate WHBB potential surface.25 In recent years, Skinner and co-workers have developed a combined electronic-structure/molecular-dynamics (ES/MD) method to simulate OH-stretch vibrational spectroscopy,26,27 and it has been successful in reproducing various OH (OD) stretch experimental spectra (e.g., IR, Raman, 2DIR, pumpprobe anisotropy, sum-frequency generation) of water in many different molecular environments (e.g., liquid water, ice Ih, air/water interface, hexamer and amorphous ices).2,26,28–43 The goal of the present work is to extend this ES/MD method to the simulation of the OH-stretch IINS of water. The rest of the paper is organized as follows. In Sec. II, we outline the general theory of IINS and the ES/MD method for IINS with small k values. In Sec. III, the calculated IINS spectra in the limit k → 0 for liquid water, supercooled water, and ice Ih are compared with experiment. In Sec. IV, the connections and differences between IR and IINS spectra are elucidated in detail. In Sec. V, we conclude.
143, 014503-1
© 2015 AIP Publishing LLC
014503-2
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
II. MIXED QUANTUM/CLASSICAL EXPRESSION FOR THE IINS SPECTRUM
For H-containing systems like water, since the IINS crosssection of hydrogen is much larger than those of other atoms including deuterium,5 it is a rather good approximation to consider only the scattering due to hydrogen in the IINS calculation. Assuming the first Born approximation and a Fermi pseudopotential between the incident neutron and the scattering atom, the IINS double-differential cross-section is proportional to the incoherent scattering function (sometimes also called the incoherent dynamical structure factor), given by44 ∞ 1 dt e−iωt I(⃗k,t), (1) S(⃗k,ω) = 2π −∞ where ⃗k is the momentum transferred during the scattering, and I(⃗k,t) is the incoherent intermediate scattering function, defined by44 1 −i k⃗ ·⃗r l (0) i k⃗ ·⃗r l (t) e , (2) I(⃗k,t) = e N l where the brackets denote an equilibrium quantum-mechanical ensemble average, N is the number of hydrogen atoms, and ⃗r l (t) is the time-dependent position operator associated with the lth hydrogen in the Heisenberg picture. Experimentally one often measures the scattering for low k (≡ |⃗k|), and defines the quantity ω2 ⃗ S( k,ω), k→ 0 k 2
G0(ω) ≡ lim
(3)
which can be called the IINS spectrum. As stated in the Introduction, it is essential to treat the OH stretch quantum-mechanically, and the ES/MD method developed by Skinner and co-workers is an accurate and efficient mixed quantum/classical approach to do so.27 In this method, the quantum-mechanical system consists of many coupled local OH stretches (also called chromophores), and its many-body vibrational ground state is denoted as |0⟩. When only the jth OH chromophore is excited to its first vibrational excited state, the vibrational state of the system is | j⟩ (note that higher excited states may be reached in nonlinear spectroscopy). The Hamiltonian (up to and including only the oneexciton band) is therefore H = H0|0⟩⟨0| + Hi j |i⟩ ⟨ j|. (4) ij
Note that in our treatment we do not include bend degrees of freedom. The system also has low-frequency modes (the bath, in this case translations and rotations of water molecules) and H0 and Hi j are operators in the bath Hilbert space. Using the completeness relation (for the ground vibrational state and the one-exciton band), 1 = |0⟩ ⟨0| + i |i⟩⟨i|, this can be rewritten as H = H0 + ~ κ i j |i⟩ ⟨ j|, (5) ij
where ~κ i j = (Hii − H0)δ i j + Hi j (1 − δ i j ).
(6)
Thus, the diagonal elements are the (typically anharmonic) local-mode transition frequencies, and the off-diagonal elements are the vibrational couplings (both intra- and intermolecular). Both frequencies and couplings are formally bath operators. H0 is the bath Hamiltonian. In what follows, the bath will become classical, and so the frequencies and couplings become functions of the bath coordinates. For a given bath configuration, these frequencies and couplings are evaluated through so-called spectroscopic “maps” based on extensive electronic-structure calculations.45 The bath coordinates evolve through a classical MD simulation, and so the frequencies and couplings fluctuate in time. In the Appendix we show that within this mixed quantum/classical approach, the quantum time-correlation function of two operators A and B can be written approximately as ⟨A0 j (0)Fj i (t)Bi0(t)⟩, (7) ⟨A(0)B(t)⟩ = ij
where now on the r.h.s. the brackets indicate a classical average, the matrix elements A0 j and Bi0 are bath variables, and the matrix F(t) satisfies the equation ˙ = iF(t)κ(t), F(t)
(8)
subject to the initial condition that Fi j (0) = δ i j . (κ i j are also bath variables.) The intermediate scattering function is therefore given approximately by 1 ⃗ ⃗ ⟨0|e−i k ·⃗r l (0)| j⟩Fj i (t)⟨i|ei k ·⃗r l (t)|0⟩ . I(⃗k,t) = N i jl
(9)
For the instantaneous position operator of hydrogen atom l, we then write ⃗l + uˆl x l , ⃗r l = R
(10)
where x l is the displacement from the equilibrium value of the H atom along the OH bond, uˆl is the unit vector of the OH ⃗l (t) is the equilibrium position of the H atom. To bond, and R first order in k, ⃗
l ⟨0|e−i k ·⃗r l | j⟩ = −ikδ jl kˆ · uˆl x 01 ,
(11)
l where x 01 is the 0-1 matrix element of the displacement coordinate for hydrogen l. For an isotropic system, G0(ω) then becomes ∞ ω2 dt e−iωt G0(ω) = 6N π −∞
l l × uˆl (0) · uˆl (t)x 01 (0)x 01 (t)Fll (t) e−|t |/2T1. (12) l
We have included the effect of the vibrational lifetime, T1, phenomenologically.27 Thus, in the end, the IINS spectrum is related to a classical time-correlation function involving the fluctuating frequencies and vibrational couplings, the fluctuating direction of each OH bond, and the fluctuating position matrix elements (which depend only weakly on the local environment).
014503-3
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
III. IINS SPECTRA FOR WATER AND ICE
To perform theoretical IINS calculations, the time evolution of the bath degrees of freedom was carried out by classical MD simulation with the third version of the E3B water model (E3B3),46 which uses the TIP4P/2005 model47 as a reference potential, and includes three-body interactions explicitly.48,49 E3B3 shows improved performance in describing some properties of water over the previous version (E3B2): for instance, better melting temperature of ice (260 K for E3B3 versus 251 K for E3B2, compared to the experimental value of 273 K), and better ice density at one atmosphere and 273 K (0.918 g/ml for E3B3 versus 0.935 g/ml for E3B2, compared to the experimental value of 0.917 g/ml).46 All simulations of liquid water, supercooled water, and ice Ih were performed in the NPT ensemble with 432-molecule simulation boxes using a modified GROMACS package 4.5.5.50–53 The pressure was maintained at 1 bar using a Parrinello-Rahman barostat54 with a coupling constant of 10.0 ps in production runs, and the temperatures were controlled by a Nosé-Hoover thermostat55,56 with a coupling constant of 5.0 ps. The geometry of the water molecule was constrained using the SETTLE algorithm,57 and the electrostatic interactions with periodic boundary conditions were treated using the particle-mesh Ewald (PME) method.58,59 The cutoff for the Lennard-Jones potential between oxygen atoms was 1.0 nm, and the longrange dispersion corrections were applied to the energy and pressure.60 A proton-disordered configuration generated by Hayward and Reimers61 was employed as the initial configuration of ice Ih. The simulation time step was 1 fs, and the trajectories were saved every 5 fs for spectral calculations. The spectroscopic maps used in the present study are from Table 1 of Ref. 45, and in order to use these TIP4P-based maps the geometry of the E3B3 (TIP4P/2005) water was adjusted to be that of TIP4P62 for the purpose of the spectral calculations. Fig. 1 displays the calculated G0(ω) (from Eq. (12)) for liquid water (top panel), supercooled water (middle panel), and (poly-crystalline) ice Ih (bottom panel), along with extrapolated results from experiment for the same systems. In the calculations, for the vibrational lifetimes we used 260 fs for ambient and supercooled water,63 and 300 fs for ice Ih.64 Note that as the melting temperature of the E3B3 model is 260 K,46 we choose 258 K for ice Ih and supercooled water in our simulations. All the theoretical and experimental spectra show essentially single peaks, except that the theoretical spectrum for ice Ih shows a shoulder on the low-frequency side of the main peak. This shoulder may also exist in experiment (it does actually exist in the 15 K experiments by Li10 but these results are not extrapolated to k = 0), but the modest frequency resolution of the experiment precludes a definitive answer. In terms of peak position, fair agreement between theory and experiment is seen for liquid and supercooled water. However, for ice Ih, the calculated spectrum is blue-shifted by about 50 cm−1 compared to the experiment. This blue shift is probably due to the fact that the frequency map overestimates the OH-stretch frequency in ice Ih by about 30-40 cm−1.45 Regarding the linewidth, our theory is somewhat too narrow for supercooled water and ice Ih, and this discrepancy might be due to the inadequacy of the experimental extrapolation, the experimental resolution, or the
FIG. 1. Experimentally extrapolated (black lines) and calculated (red lines) G0(ω) for OH stretch in liquid water (top panel), supercooled water (middle panel), and ice Ih (bottom panel). Reference numbers for experiments are: (a) Ref. 16, (b) Ref. 11, (c) Ref. 14. All the spectra are normalized to have the same height of 1.
fact that our simulation temperatures (258 K) are lower than corresponding experimental ones (∼270 K). Overall, however, the agreement between theory and experiment is reasonable, especially in reproducing the red peak-shifts from liquid water, to supercooled water, to ice Ih (from top to bottom). IV. COMPARISON OF IINS AND IR SPECTRA
OH-stretch IR and Raman spectroscopies have been used to study liquid water and ice for years, and there are a number of experimental spectra available in the literature (see citations in Refs. 2 and 34). From the top and bottom panels of Fig. 1, one might already realize one difference between experimental IR line-shapes and G0(ω) for water: the IR spectra are peaked at lower frequencies (about 3400 cm−1 for liquid water at room temperature,1 and about 3250 cm−1, for ice at about 270 K34,65) compared to the corresponding G0(ω) (about 3450 cm−1 for liquid water,12,16 and about 3350 cm−1 for ice14,16). In this section, we will review the theoretical results for IR spectra, within the same mixed quantum/classical approach, and then highlight the two major differences between IR and IINS spectra. The IR spectrum for an isotropic system is related to the Fourier transform of the dipole-dipole time-correlation function,44 ∞ 1 ⃗ ⃗ I(ω) = dt e−iωt M(0) · M(t) , (13) 2π −∞ ⃗ is the total dipole of the sample (in the ground elecwhere M tronic state). The total dipole can be expanded to first order in the local-mode displacements. Within the bond-dipole approximation the required matrix element (according to Eq. (7)) is therefore ⃗ | j⟩ = µ′j uˆ j x j , ⟨0| M 01
(14)
014503-4
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
where µ′ is the dipole derivative, and the other symbols were defined previously. The spectrum, then, becomes2,3 ∞ 1 j i I(ω) = µ′j (0)µ′i (t)uˆ j (0) · uˆi (t)x 01 (0)x 01 (t)Fj i (t) e−|t |/2T1. dt e−iωt 2π −∞ ij
(15)
This expression involves a double sum over chromophores, and so there are two distinct contributions, when i = j and when i , j. Explicitly, they are ∞ 1 j j −iωt ′ ′ Iinc(ω) = dt e µ j (0)µ j (t)uˆ j (0) · uˆ j (t)x 01(0)x 01(t)Fj j (t) e−|t |/2T1 (16) 2π −∞ j and 1 Icoh(ω) = 2π
∞
dt e −∞
−iωt
µ′j (0)µ′i (t)uˆ j (0)
·
j i uˆi (t)x 01 (0)x 01 (t)Fj i (t)
e−|t |/2T1.
(17)
i, j
The first term is closely related to Eq. (12) for IINS, and so we call it “incoherent.” We call the second term “coherent” for two reasons: (1) it is similar to the major contribution to coherent neutron scattering; (2) this decomposition has been used previously to analyze IR line shapes,3 and this second term is responsible for “coherent” vibrational energy transfer. In Fig. 2, the calculated IR spectra (solid black lines) and their incoherent and coherent contributions (dashed and dotted black lines, respectively) for supercooled water (right panel) and ice Ih (left panel) at 258 K are shown in the top panels. One sees that the full line shape arises from a significant cancellation of the two terms. The coherent term would be zero in the absence of vibrational coupling and energy transfer.3 In the same figure, bottom panels, we show the calculated corresponding G0(ω) (red solid lines). It is evident that in each case the incoherent IR line shape Iinc(ω) closely resembles G0(ω). In fact, Eq. (12) for IINS and Eq. (16) for the incoherent contribution to IR are the same except for the factors of the dipole derivative, µ′, in the latter. µ′ is known to depend sensitively on the molecular environment, and stronger hydrogen bonding leads to a larger value.66 This “non-Condon” effect
has been shown to be important for various kinds of IR spectroscopy in liquid water. Since stronger hydrogen bonding corresponds to lower frequency, non-Condon effects tend to shift the spectrum to the red. Indeed, in Fig. 2 the incoherent contribution to the IR spectrum is red-shifted (slightly, by about 20 cm−1) compared to the IINS spectrum for supercooled water. In ice, however, all molecules have roughly the same hydrogen-bonding environment (they all have four hydrogen bonds), and so the dipole derivative is more or less constant. In this case, G0(ω) ∼ ω2 Iinc(ω), and indeed, in Fig. 2 we see that is the case (the very weak ω2 dependence can be ignored). To illustrate the relation between the VDOS and G0(ω), we consider the system in the inhomogeneous limit, when bath dynamics becomes unimportant.2,67 We also take T1 → ∞. In this limit, each of the fluctuating quantities at time t can be replaced by its value at time 0, and so, for example, we have ∞ l G0(ω)/ω2 ∼ dt e−iωt (0)2[eiκ(0)t ]ll . (18) x 01 −∞
l
FIG. 2. Calculated OH stretch IR and IINS spectra for ice Ih (left panel) and supercooled water (right panel) at 258 K. IR line shapes (solid black lines) with incoherent (dashed black lines, from Eq. (16)) and coherent (dotted black lines, from Eq. (17)) contributions are shown in the top, and IINS spectra G0(ω) (red lines, from Eq. (12)) are in the bottom. Frequency distributions (solid blue lines), calculated from Eq. (21), are also shown in the bottom panel.
014503-5
L. Shi and J. L. Skinner
J. Chem. Phys. 143, 014503 (2015)
The matrix elements fluctuate only weakly, and therefore to an excellent approximation can be taken outside the average, in which case, ∞ 2 G0(ω)/ω ∼ dt e−iωt Tr[eiκ(0)t ] . (19) −∞
Since this result is now independent of the time origin, the argument of κ can be dropped. Further, κ can be diagonalized by an orthogonal matrix, QT κQ = λ,
(20)
where λ is diagonal, and the eigenvalues λ l are the frequencies of the vibrational eigenmodes. Therefore, it is easy to see that in this limit the IINS spectrum is 2 G0(ω)/ω ∼ δ(ω − λ l ) ≡ f (ω). (21) l
f (ω) is the VDOS. Note that this is not simply a harmonic result; rather, this is the VDOS even for arbitrarily large localmode anharmonicity. The calculated f (ω) for supercooled water and ice Ih at 258 K are shown as solid blue lines in Fig. 2, and the similarity between G0(ω) and f (ω) demonstrates again that IINS is useful in extracting the VDOS, although G0(ω) is sharpened from f (ω) due to motional narrowing. In contrast, the total IR line shape (solid black lines) is quite different from the VDOS both in terms of peak position and linewidth due to non-Condon effects and the coherent contribution.
ACKNOWLEDGMENTS
This work was supported by NSF Grant No. CHE1058752.
APPENDIX: MIXED QUANTUM/CLASSICAL EXPRESSION FOR QUANTUM TIME-CORRELATION FUNCTIONS
The equilibrium quantum time-correlation function for two arbitrary operators A and B can be written as ⟨A(0)B(t)⟩ = Tr[e−β H Aei H t /~ Be−i H t /~]/Tr[e−β H ],
(A1)
where in this case the Hamiltonian is given by Eq. (5) (and β = 1/k BT). The trace is over the bath states, and over the ground and one-exciton states. For these high-frequency vibrations ~κ ii ≫ k BT, and so the vibrational part of the trace involves only the ground state. Moreover, for the frequencies (in the eventual Fourier transform) of interest ω ∼ κ ii , this becomes A0 j ⟨ j|ei H t /~|i⟩Bi0e−i H0t /~]/ ⟨A(0)B(t)⟩ = Tr b [e−β H0 ij
Tr b [e
−β H0
],
(A2)
where Tr b indicates a trace over bath states. Bi0 is a bath operator, and so its time dependence is Bi0(t) = ei H0t /~ Bi0e−i H0t /~. With this definition, this becomes ⟨A(0)B(t)⟩ = Tr b [e−β H0 A0 j Fj i (t)Bi0(t)]/Tr b [e−β H0], ij
(A3) where
V. CONCLUDING REMARKS
In recent years, due to the availability of high-energy neutron sources, OH-stretch IINS spectra for water in condensed phases have been measured, and OH-stretch VDOS have been extracted from these experimental measurements, providing a good basis to evaluate modeling methods of OHstretch vibrational spectroscopy. In this work, we extended the ES/MD method, developed by Skinner and co-workers for optical vibrational spectroscopy (e.g., IR and Raman), to simulating OH-stretch IINS spectra at small k values. The reasonable agreement between theory and experiment for liquid water, supercooled water, and ice Ih further validates the ES/MD method, which has been quite successful in modeling various optical vibrational spectroscopies. We then go on to analyze the differences between IINS and IR spectra. We find there are two: first, the IR spectrum involves the dipole derivative, which tends to weight lower frequencies more; second and more importantly, the IR spectrum involves two terms—one that is analogous to incoherent neutron scattering and another that is analogous to coherent scattering. This second term involves the correlations and couplings between different molecules, and is signed (can be positive or negative). Adding this term to the “incoherent” contribution can dramatically change the spectrum. For both of these reasons the IR spectrum is more difficult to interpret than the IINS spectrum, which to a reasonably good approximation is just the VDOS.
Fj i (t) = ⟨ j|ei H t /~e−i H0t /~|i⟩.
(A4)
From Eq. (5), we have F(t) = ei(H0+~κ)t /~e−i H0t /~.
(A5)
From this it is clear that ˙ = iF(t)κ(t), F(t)
(A6)
where κ(t) = ei H0t /~ κe−i H0t /~. At this point we take the classical limit, such that F(t) and κ(t) are matrices in the oneexciton Hilbert space, but the matrix elements themselves are classical variables. In this limit the trace over the bath, weighted by the bath density operator, can be replaced by the classical equilibrium statistical mechanical average, recovering Eq. (7). 1J.
E. Bertie and Z. Lan, Appl. Spectrosc. 50, 1047 (1996). M. Auer and J. L. Skinner, J. Chem. Phys. 128, 224511 (2008). 3M. Yang and J. L. Skinner, Phys. Chem. Chem. Phys. 12, 982 (2010). 4P. C. H. Mitchell, S. F. Parker, A. J. Ramirez-Cuesta, and J. Tomkinson, Vibrational Spectroscopy with Neutrons With Applications in Chemistry, Biology, Materials Science and Catalysis (World Scientific Publishing Co. Pte. Ltd., Singapore, 2005). 5S. F. Parker, D. Lennon, and P. W. Albers, Appl. Spectrosc. 65, 1325 (2011). 6C. M. Brown, A. J. Ramirez-Cuesta, M. R. Johnson, and V. García-Sakai, Chem. Phys. 427, 1 (2013). 7S.-H. Chen, K. Toukan, C.-K. Loong, and D. L. Price, Phys. Rev. Lett. 53, 1360 (1984). 8C. Andreani, P. Bosi, F. Sacchetti, and C. K. Loong, J. Chem. Phys. 83, 750 (1985). 2B.
014503-6 9K.
L. Shi and J. L. Skinner
Toukan, M. A. Ricci, S.-H. Chen, C.-K. Loong, D. L. Price, and J. Teixeira, Phys. Rev. A 37, 2580 (1988). 10J.-C. Li, J. Chem. Phys. 105, 6733 (1996). 11M. A. Ricci, M. Nardone, A. Fontana, C. Andreani, and W. Hahn, J. Chem. Phys. 108, 450 (1998). 12D. D. Klug, C. A. Tulk, E. C. Svensson, and C.-K. Loong, Phys. Rev. Lett. 83, 2584 (1999). 13Y. Wang, A. I. Kolesnikov, S. L. Dong, and J. C. Li, Can. J. Phys. 81, 401 (2003). 14E. Perelli-Cippo, G. Gorini, M. Tardocchi, C. Andreani, A. Pietropaolo, R. Senesi, N. J. Rhodes, and E. M. Schooneveld, Appl. Phys. A 83, 453 (2006). 15C. Andreani, G. Romanelli, and R. Senesi, Chem. Phys. 427, 106 (2013). 16R. Senesi, D. Flammini, A. I. Kolesnikov, E. D. Murray, G. Galli, and C. Andreani, J. Chem. Phys. 139, 074504 (2013). 17E. J. Heller, J. R. Reimers, and G. Drolshagen, Phys. Rev. A 36, 2613 (1987). 18T. D. Hone and G. A. Voth, J. Chem. Phys. 121, 6412 (2004). 19I. R. Craig and D. E. Manolopoulos, Chem. Phys. 322, 236 (2006). 20J. Liu and W. H. Miller, J. Chem. Phys. 128, 144511 (2008). 21S. Bratos, G. Tarjus, M. Diraison, and J.-C. Leicknam, Phys. Rev. A 44, 2745 (1991). 22S. Bratos, M. Diraison, G. Tarjus, and J.-C. Leicknam, Phys. Rev. A 45, 5556 (1992). 23G. Tarjus, M. Diraison, J.-C. Leicknam, and S. Bratos, J. Mol. Struct. 296, 199 (1993). 24M. Diraison, J.-C. Leicknam, G. Tarjus, and S. Bratos, Phys. Rev. E 50, 2689 (1994). 25H. Liu, Y. Wang, and J. M. Bowman, J. Phys. Chem. B 117, 10046 (2013). 26S. A. Corcelli, C. P. Lawrence, and J. L. Skinner, J. Chem. Phys. 120, 8107 (2004). 27J. L. Skinner, B. M. Auer, and Y.-S. Lin, Adv. Chem. Phys. 142, 59 (2009). 28S. A. Corcelli and J. L. Skinner, J. Phys. Chem. A 109, 6154 (2005). 29B. M. Auer, R. Kumar, J. R. Schmidt, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 104, 14215 (2007). 30Y. Ni and J. L. Skinner, J. Chem. Phys. 141, 024509 (2014). 31F. Li and J. L. Skinner, J. Chem. Phys. 132, 204505 (2010). 32F. Li and J. L. Skinner, J. Chem. Phys. 133, 244504 (2010). 33F. Li and J. L. Skinner, J. Chem. Phys. 134, 099901 (2011). 34L. Shi, S. M. Gruenbaum, and J. L. Skinner, J. Phys. Chem. B 116, 13821 (2012). 35L. Shi and J. L. Skinner, J. Phys. Chem. B 117, 15536 (2013). 36L. Shi, F. Li, and J. L. Skinner, J. Chem. Phys. 140, 244503 (2014). 37I. V. Stiopkin, C. Weeraman, P. A. Pieniazek, F. Y. Shalhout, J. L. Skinner, and A. V. Benderskii, Nature 474, 192 (2011). 38P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Am. Chem. Soc. 133, 10360 (2011).
J. Chem. Phys. 143, 014503 (2015) 39P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Chem. Phys. 135, 044701
(2011). Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 110, 1992 (2013). 41C. J. Tainter, Y. Ni, L. Shi, and J. L. Skinner, J. Phys. Chem. Lett. 4, 12 (2013). 42C. J. Tainter, L. Shi, and J. L. Skinner, J. Chem. Phys. 140, 134503 (2014). 43T. L. C. Jansen, B. M. Auer, M. Yang, and J. L. Skinner, J. Chem. Phys. 132, 224503 (2010). 44D. A. McQuarrie, Statistical Mechanics (Harper and Row, New York, 1976). 45S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, J. Chem. Theory Comput. 9, 3109 (2013). 46C. J. Tainter, L. Shi, and J. L. Skinner, J. Chem. Theory Comput. 11, 2268 (2015). 47J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005). 48R. Kumar and J. L. Skinner, J. Phys. Chem. B 112, 8311 (2008). 49C. J. Tainter, P. A. Pieniazek, Y.-S. Lin, and J. L. Skinner, J. Chem. Phys. 134, 184501 (2011). 50B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). 51D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen, J. Comput. Chem. 26, 1701 (2005). 52E. Lindahl, B. Hess, and D. van der Spoel, J. Mol. Model. 7, 306 (2001). 53H. J. C. Berendsen, D. van der Spoel, and R. van Drunen, Comput. Phys. Commun. 91, 43 (1995). 54M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982). 55S. Nosé, Mol. Phys. 52, 255 (1984). 56W. G. Hoover, Phys. Rev. A 31, 1695 (1985). 57S. Miyamoto and P. A. Kollman, J. Comput. Chem. 13, 952 (1992). 58T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993). 59U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995). 60M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford, 1987). 61J. Hayward and J. Reimers, J. Chem. Phys. 106, 1518 (1997). 62W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 63A. J. Lock and H. J. Bakker, J. Chem. Phys. 117, 1708 (2002). 64R. L. A. Timmer and H. J. Bakker, J. Phys. Chem. A 114, 4148 (2010). 65S. G. Warren and R. E. Brandt, J. Geophys. Res. 113, D14220 doi:10.1029/2007jd009744 (2008). 66J. R. Schmidt, S. A. Corcelli, and J. L. Skinner, J. Chem. Phys. 123, 044513 (2005). 67B. M. Auer and J. L. Skinner, J. Chem. Phys. 127, 104105 (2007). 40Y.
