IR and SFG vibrational spectroscopy of the water bend in the bulk liquid and at the liquid-vapor interface, respectively Yicun Ni and J. L. Skinner

Citation: J. Chem. Phys. 143, 014502 (2015); doi: 10.1063/1.4923462 View online: http://dx.doi.org/10.1063/1.4923462 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 Communication: Vibrational sum-frequency spectrum of the air-water interface, revisited J. Chem. Phys. 145, 031103031103 (2016); 10.1063/1.4958967 IR spectra of water droplets in no man’s land and the location of the liquid-liquid critical point J. Chem. Phys. 145, 124509124509 (2016); 10.1063/1.4963736 Ultrafast pump-probe and 2DIR anisotropy and temperature-dependent dynamics of liquid water within the E3B model J. Chem. Phys. 141, 024509024509 (2014); 10.1063/1.4886427 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 J. Chem. Phys. 144, 214501214501 (2016); 10.1063/1.4952991

THE JOURNAL OF CHEMICAL PHYSICS 143, 014502 (2015)

IR and SFG vibrational spectroscopy of the water bend in the bulk liquid and at the liquid-vapor interface, respectively Yicun Ni and J. L. Skinner Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, USA

(Received 22 May 2015; accepted 23 June 2015; published online 7 July 2015) Vibrational spectroscopy of the water bending mode has been investigated experimentally to study the structure of water in condensed phases. In the present work, we calculate the theoretical infrared (IR) and sum-frequency generation (SFG) spectra of the HOH bend in liquid water and at the water liquid/vapor interface using a mixed quantum/classical approach. Classical molecular dynamics simulation is performed by using a recently developed water model that explicitly includes three-body interactions and yields a better description of the water surface. Ab-initio-based transition frequency, dipole, polarizability, and intermolecular coupling maps are developed for the spectral calculations. The calculated IR and SFG spectra show good agreement with the experimental measurements. In the theoretical imaginary part of the SFG susceptibility for the water liquid/vapor interface, we find two features: a negative band centered at 1615 cm−1 and a positive band centered at 1670 cm−1. We analyze this spectrum in terms of the contributions from molecules in different hydrogen-bond classes to the SFG spectral density and also compare to SFG results for the OH stretch. SFG of the water bending mode provides a complementary picture of the heterogeneous hydrogen-bond configurations at the water surface. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4923462]

I. INTRODUCTION

Vibrational spectroscopy is an excellent experimental method to study the microscopic structure and dynamics of condensed-phase water. In the mid-infrared (IR) frequency region, most work has focused on the intramolecular OHstretch motion because this vibrational mode is very sensitive to its local hydrogen-bonding environment.1 This technique can be particularly illuminating if dilute HOD in either H2O or D2O is studied, because in this case, the complicating effects of intramolecular and intermolecular vibrational couplings can be greatly reduced.1–10 Vibrational spectroscopy can also be used to study properties at interfaces. Sum-frequency-generation (SFG) spectroscopy, a surface-sensitive technique, has been applied to study the structure of water at interfaces.11–22 Tian and Shen first applied heterodyne-detected SFG spectroscopy to the dilute HOD in D2O system to study the hydrogen-bond structure of the water liquid/vapor interface.23 The measured spectrum shows three major features: a sharp positive peak at about 3700 cm−1 corresponding to the upward-pointing dangling OH bonds and negative and positive bands at 3500 cm−1 and 3300 cm−1, respectively, arising from hydrogen-bonded OH groups. The latter two peaks were interpreted, respectively, as “water-like” molecules with downward-pointing OH bonds and “ice-like” molecules with upward-pointing OH bonds by these authors.23,24 Others have formulated alternative interpretations.25–30 In a recent letter,31 we have attempted to provide a unified interpretation of the OH-stretch spectroscopy of HOD in different D2O environments: the hexamer cage, ice, the bulk liquid, and the liquid/vapor interface. The last two are particularly relevant here. In our analysis, we needed to consider 0021-9606/2015/143(1)/014502/12/$30.00

different classes of water molecules, which we labeled by the total number of hydrogen bonds from a molecule and by the number of “donor” hydrogen bonds.32 For example, a 3 D molecule is one with three hydrogen bonds, two of which are donors (D stands for “double,” S for “single,” and N for “nodonors”). We then decomposed the appropriate spectral densities in terms of hydrogen-bonded donor-acceptor pairs, where the HOD molecule donates its H atom for hydrogen bonding with a D2O molecule, and also in terms of HOD molecules without a hydrogen bond to the H (“free-OH molecules”). We found that for the IR spectrum of the bulk liquid, the main contributing donor-acceptor pairs were 4 D -4 D , 4 D -3 D , and 3 D -4 D , and there were also some 2S and 3S free-OH molecules. All of these contributions are very broad and overlapping, leading to an IR spectrum with a single peak. In terms of the SFG spectrum, we found that the free-OH groups on 2S and 3S molecules led to the sharp positive peak at 3700 cm−1, the 2S -4 D pair was the main contributor to the negative peak at 3500 cm−1, and the 4 D -2S pair was the main contributor to the positive peak at 3300 cm−1. Many other donor-acceptor pairs are present at the surface, but they produce smaller contributions to the SFG spectrum, often due to cancellation. Although the bend vibration of water is less sensitive to molecular environment than the stretch, it nonetheless provides a useful, and complementary, probe of water structure. In the gas phase, the frequency of the water bending mode is 1595 cm−1.33 It is blue-shifted to a higher frequency in either the bulk liquid (1645 cm−1 34–38) or in ice (1670 cm−1 39). Vibrational spectroscopy has not been performed as often for the water bend as it has for the OH stretch. Some of the reasons are the smaller transition dipole moment of the bend and its weaker frequency dependence on its local environment. However, the

143, 014502-1

© 2015 AIP Publishing LLC

014502-2

Y. Ni and J. L. Skinner

bending mode has its own advantages: first, there is only one bending mode per molecule, which means there is essentially no intramolecular vibrational coupling (for the fundamental). Second, due to its smaller transition dipole moment, the bend line shapes should not be affected too much by intermolecular vibrational coupling.40 In addition, it is important to study the bending mode since it is a major channel for vibrational energy relaxation processes from higher frequency intramolecular modes, i.e., OH stretch, to intermolecular vibrational modes such as the libration and hydrogen-bond stretch and bend.41–52 Finally, studying the bending mode can help researchers understand the nature of the Fermi resonance in water.5,53–56 Experimentally, Falk measured the IR spectrum of the water bend in different liquids and solids.33 The temperature dependence of the IR spectrum of the bend mode in liquid water has also been studied; a slight blue-shift has been observed as temperature decreases.41,57 This blue-shift is interpreted as being caused by a stronger hydrogen-bond network between water molecules at lower temperature. Pavlovic et al. measured Raman spectra of the bend in liquid water with different polarizations.38 By decomposing the line shape, these authors found that two components were needed to have a satisfactory fit to the isotropic spectrum, which in their view supported the existence of two structures in liquid water, corresponding to strongly and weakly hydrogen-bonded molecules. Vinaykin and Benderskii40 and Nagata et al.58 have recently measured homodyne-detected SFG spectra of the water liquid/vapor interface, in which two features, one at 1620 cm−1 and the other one at 1685 cm−1, are found. These have been assigned to molecules with a free OH and to molecules with two donor hydrogen bonds, respectively.40,58 In terms of theoretical vibrational spectroscopy of liquid water, the simplest approach uses a classical molecular dynamics (MD) simulation with flexible water molecules and typically uses a harmonic quantum-correction factor to account for quantum effects due to the high frequencies of some of the modes (as a result of the light H atoms).58–68 Such an approach can yield IR or SFG spectra over the entire frequency range from 0-4000 cm−1, encompassing the low-frequency modes, librations, the bend, bend-libration combinations, and the OHstretch region. Specifically, results for the SFG spectrum in the bend region have recently appeared.58,61 One can also account for quantum effects more accurately and directly with path integrals, centroid MD, or ring-polymer MD.69–72 As another alternative, one can perform ab-initio MD simulations, where the Born-Oppenheimer surface is computed on the fly.9,73–78 Quantum effects can be incorporated into these calculations as well. While some of these methods have led to quite good results, in others the vibrational frequencies are not so accurate, either because vibrational anharmonicities are not taken into account properly, zero-point energies are neglected, or because of inaccuracies in the forms of the intramolecular potentials chosen in the empirical models or in the density functional theory (DFT) used in the ab-initio MD calculations. Focusing only on the OH-stretch region, we have developed an alternative mixed quantum/classical approach, where the low-frequency modes are treated classically through MD simulation, and the OH stretches are treated quantummechanically through a vibrational exciton Hamiltonian in the

J. Chem. Phys. 143, 014502 (2015)

local-mode basis.79,80 The diagonal elements of the exciton Hamiltonian are the local-mode 0-1 anharmonic transition frequencies, while the off-diagonal elements are the intramolecular and intermolecular couplings between local modes. These exciton matrix elements are all obtained from electronic-structure (ES)-based “maps.”81 In addition, one needs the local-mode transition dipole moments, which are obtained from an additional map.81 With this approach, we have calculated a variety of different vibrational spectroscopies (IR, Raman, 2DIR, SFG, and 2DSFG) for liquid water, ice, amorphous ices, the liquid/vapor and other interfaces, the water hexamer, salt solutions, and confined water.25,31,79,82–92 The purpose of this paper is to develop a similar mixed quantum/classical method for the water bend, with applications to IR and SFG spectroscopy of the bulk liquid and liquid/vapor interface, respectively. Our goal is to develop an approach that can reproduce experiment and then to use our simulation results to provide information about molecular structure in the bulk liquid and at the water liquid/vapor interface. Such information is complementary to that obtained from OH-stretch spectroscopy, since the bend and local-mode stretch have different molecule-fixed transition dipole moments, and the transition frequencies are sensitive in different ways to local hydrogen bonding. Taken together, bend and stretch spectroscopies can provide a more complete molecular picture of structure in the bulk and at the surface. The paper is organized as follows. In Section II, the theoretical background for spectral calculations is described. In Section III, we discuss new spectroscopic maps for the bend mode. In Section IV, we present our calculated IR and SFG spectra of the water bend and make comparisons to those from experiments. Section V analyzes the IR and SFG spectra by decomposing the spectral densities into contributions from different hydrogen-bond classes. We also make the connection to previous results for OH-stretch spectroscopy. We conclude in Section VI.

II. CALCULATION OF LINE SHAPES A. Mixed quantum/classical formalism

Both IR and SFG line shapes can be written in terms of the Fourier-Laplace transform of quantum time-correlation functions (TCFs). For example, the IR absorption line shape is expressed as93  I p (ω) ∝ Re

0

∞

dt eiωt ⟨ pˆ · ⃗µ(t) ⃗µ(0) · p⟩ ˆ ,

(1)

where pˆ is the unit vector of the polarization of the excitation light. ⃗µ(t) is the time-dependent dipole operator (in the ground electronic state) of the system at time t, and the brackets indicate a quantum equilibrium ensemble average. Within the mixed quantum/classical approach, for a system of coupled chromophores, the IR line shape can be written as79,94  ∞ 

 I p (ω) ∝ Re dt eiωt mi p (t)Fi j (t)m j p (0) e−t /2T1, 0

ij

(2)

014502-3

Y. Ni and J. L. Skinner

J. Chem. Phys. 143, 014502 (2015)

where mi p (t) is the component of the transition dipole for localmode chromophore i along the p direction at time t. Fi j (t) are the elements of the propagator matrix F(t), which satisfies the equation of motion, ˙ = −iF(t)κ(t), F(t)

(3)

subject to the initial condition that Fi j (0) = δ i j and with κ i j (t) = ωi (t)δ i j + ωi j (t)(1 − δ i j ),

(4)

where κ(t) is the exciton Hamiltonian matrix (divided by ~) whose diagonal elements are the fluctuating transition frequencies, ωi (t), and whose off-diagonal elements are the fluctuating couplings between different chromophores, ωi j (t). The angular brackets now denote a classical equilibrium ensemble average. As before, the effect of the vibrational lifetime, T1, is included phenomenologically as shown; the value of T1 will be taken from experiment. If the system contains only a single isolated chromophore, Eq. (2) can be simplified as95  ∞ I p (ω) ∝ Re dt eiωt 0  (  t ) × m p (t)m p (0) exp −i dτω(τ) e−t /2T1. (5) 0

Especially in this case, it is easier and often more instructive to consider an inhomogeneous approximation to Eq. (5). If the dynamics of the system is sufficiently slow, then ω(τ) can be replaced by ω(0) and m p (t) can be replaced by m p (0). If, in addition, we neglect lifetime effects (by setting T1 = ∞), then

 I p (ω) ∝ m p (0)2δ(ω − ω(0)) , (6) which is the “spectral density.”26,96 Within the Condon approximation (the magnitude of the transition dipole is constant, and so m p (0) is uncorrelated with ω(0)), this is simply proportional to the distribution of frequencies P(ω) = ⟨δ(ω − ω(0))⟩. Similarly, in a heterodyne-detected SFG experiment, the imaginary part of the second-order susceptibility is measured. Within the mixed quantum/classical approach, for a system with coupled chromophores, the tensor components of the (complex) resonant SFG susceptibility are26,97  ∞  R χ pqr (ω) ∝ i dt eiωt ⟨ai pq (t)Fi j (t)m jr (0)⟩e−t /2T1, 0

ij

(7) where ai pq (t) is the pq component (in the lab frame) at time t of the transition polarizability tensor for chromophore i. And for an isolated chromophore, the above equation can be simplified as26,96  ∞ χ Rpqr (ω) ∝ i dt eiωt 0  (  t ) × a pq (t)mr (0) exp −i dτω(τ) e−t /2T1. (8) 0

26,96

The SFG spectral density is

 Im( χ Rpqr (ω)) ∝ a pq (0)mr (0)δ(ω − ω(0)) .

(9)

III. SPECTROSCOPIC MAPS A. ES/MD approach and the frequency map

In the above formulas, one sees that IR and SFG spectra of the water bend can be calculated if the trajectories of the fluctuating transition frequencies and couplings, and transition dipole moments and polarizabilities are known. To do this, we follow the ES/MD approach proposed by Skinner and coworkers to calculate the spectra of the OH-stretch region.81,94,95 The philosophy of this method is to build empirical relationships, called “maps,” between the fluctuating quantities and some collective variable of the system. In the case of the OH stretch, such a variable is the projection of the electric field (due to the point charges on the other water molecules) on the OH bond of interest. Therefore, to calculate the vibrational spectra of the water bend, similar maps are needed. To this end, following previous studies,81,95 1500 water clusters are generated from a MD simulation of TIP4P water at 300 K. To construct a cluster, one first selects a water molecule at random as the central water molecule and defines its HOH angle as the angle of interest. Any surrounding molecule having its oxygen atom within 4 Å of either hydrogen atoms on the central molecule is explicitly included in the cluster, and mostly for historical reasons, any molecules having its oxygen within 7.831 Å (but larger than 4 Å) of either one of those hydrogen atoms is represented by its TIP4P point charges. To keep the configurational independence among those clusters, at most one cluster is extracted from each snapshot and this is repeated every picosecond in the process of the MD trajectory until the required number of clusters is collected. Electronicstructure calculations will be performed on these clusters (see below). The full vibrational Hamiltonian for a single molecule can be expressed as98 Hˆ = Tˆ + Vˆ ~2 ∂ 2 ~2 ∂2 ~2 ∂ 2 − − z =− 2µOH ∂r 12 2µOH ∂r 22 mO ∂r 1∂r 2 ~2  1 ∂ ∂  1 2z − + − , (1 − z 2) 2 2 4 µOH r 1 µOH r 2 mOr 1r 2 ∂z ∂z + ~2 1 ∂ 1 ∂
∂ − + (1 − z 2)1/2 (1 − z 2)1/2 mO r 1 ∂r 2 r 2 ∂r 1 ∂z ˆ + V (r 1,r 2, z), (10) where r 1 and r 2 are the two OH distances, respectively, z ≡ cos θ and θ is the HOH angle, mO is the mass of an oxygen atom, µOH is the reduced mass for oxygen and hydrogen atoms, and [, ]+ is the anti-commutator. Considering a water molecule brought from gas to liquid phases, its OH bondlength increases because of the formation of hydrogen bonds. Therefore, although we set the OH bond to be rigid in our bend-mode frequency calculation, we use an elongated OH bond length instead of its gas-phase value. To calculate the OH bond-length in the liquid phase, we perform constrained geometrical optimization on the central molecule in each of the above 1500 clusters, where the configurations of surrounding molecules, and the center of mass and the orientation of the central molecule are kept unchanged. The optimized bond

014502-4

Y. Ni and J. L. Skinner

length, which is 0.9875 Å on average, is used in the HOH bend Hamiltonian. Under this (rigid-stretch) approximation, the HOH bend Hamiltonian can be expressed as 2z ∂ ∂  ~2  2 − , (1 − z 2) + Vˆ (z), Hˆ = Tˆ + Vˆ = − 4 µOH r 2 mOr 2 ∂z ∂z + (11) where r is the above OH bond-length and other symbols have the same meanings as in Eq. (10). To scan the (onedimensional) potential energy surface (PES) for the HOH angle, Vˆ (z), in each water cluster the HOH angle of the central water molecule is stepped from 60◦ to 140◦ with a grid spacing of 2◦. The center of mass and the orientation of the central water molecule are again kept fixed during the variation of the angle. Single-point energy electronic-structure calculations are performed by the Gaussian09 software package at the theory level of B3LYP/6−311++G∗∗. From the PES and Eq. (11), the 1-0 bending frequency ω10 was then computed using the discrete variable representation (DVR) method of Colbert and Miller.99 When this procedure is applied to an isolated gas-phase water molecule, it yields a HOH bend frequency of 1564 cm−1, slightly below the experimental number of 1595 cm−1. To compensate for this, all calculated frequencies are scaled by a factor of 1.0198. These calculated frequencies are shown as dots in Fig. 1, panel (a). For completeness (although we do not need these results in this paper), we show the anharmonicities (∆ = ω10 − ω21) in panel (b). As for the OH-stretch mode, we need to correlate these frequencies with a collective coordinate of the system. In the case of the HOH bend, it is natural to choose the sum of the

FIG. 1. DFT HOH bend-fundamental frequency ω 10 (panel (a)) and anharmonicity ∆ (panel (b)) versus the electric field E (blue dots). The dashed black lines show the map results, and the black dots are the DFT results for the water monomer.

J. Chem. Phys. 143, 014502 (2015)

components of the electric fields perpendicular to the two OH bonds and pointing towards the HOH bisector, evaluated at the H atoms, from the TIP4P point-charge representation of all the molecules in the cluster surrounding the central H2O molecule. In Fig. 1, the fundamental frequencies and anharmonic shifts are plotted against this electric-field component sum, and one finds a rough linear dependence. To determine the frequency map, we perform a least-squares fit and the results are shown as the dashed line in Fig. 1. The root-mean squared deviation (RMSD) of the points from the fit is 23 cm−1. In an effort to reduce this quantity, we did try other more sophisticated collective coordinates, but they did not lead to significant improvement. B. Transition dipole and polarizability

The projection of the HOH transition dipole moment on the laboratory frame pˆ axis can be expressed as95 m p = pˆ · ⃗µ10 ≈ pˆ · ⃗µ′θ 10 ≈ µ′θ 10 uˆ · p, ˆ

(12)

where µ′ is the magnitude of the dipole derivative, θ 10 is the 1-0 transition matrix element of the HOH bend coordinate, and uˆ is the unit vector in the direction of the molecule-fixed transition dipole, in this case opposite to the HOH bisector. The transition matrix element θ 10 (and for completeness θ 21) can be readily calculated from the DVR scheme, and they depend on the bath coordinates just as the transition frequency does. In Fig. 2, we show linear fits of the calculated matrix elements to the calculated frequencies. The magnitude of the dipole derivative can also be obtained from the electronic-structure calculations

FIG. 2. Calculated matrix element θ 10 vs ω 10 (panel (a)) and θ 21 vs ω 21 (panel (b)).

014502-5

Y. Ni and J. L. Skinner

J. Chem. Phys. 143, 014502 (2015)

for each cluster by using numerical finite differences for the total dipole moment in the vicinity of the bottom of the PES. We find that µ′ for the bend mode is not very sensitive to the environment, as it is for the OH stretch mode (data not shown), which means that the Condon approximation is applicable for the bend. Thus, we use the average value of the magnitude of the bend dipole derivative in all clusters, which is 0.321 63 a.u. (atomic units). The weak non-Condon effect of the water bend was also remarked previously.40 The transition polarizability of the bend is taken to be its gas-phase value, and so it is also independent of the local environment. We write a pq = α ′pq θ 10.

(13)

In the molecular frame of a water molecule, we set the direction of the HOH bisector as the x axis, the direction perpendicular to the HOH plane as the y axis, and the H-H vector as the z axis. For this molecular frame, the polarizability derivative tensor is diagonal. It is well known that the calculation of polarizabilities requires very large basis sets.100 We use B3LYP/daug-ccpVTZ to calculate the polarizability derivative tensor of the gas-phase molecule and obtain α ′x x = −0.220 589 a.u., α ′y y = 0.544 310 a.u., and α z′ z = 1.461 042 a.u., respectively. C. Intermolecular coupling

For the vibrational coupling between HOH bend chromophores on different molecules, the simplest model uses the transition-dipole approximation.81,94 In this model, the interaction between two chromophores originates from dipole-dipole interactions and the position of the transition dipole moment on the chromophore can be obtained by fitting to electronicstructure calculations. This methodology has been successfully applied to the calculation of the OH-stretch vibrational coupling.81,94 However, when applying this model to the bend mode, we find, surprisingly, that the calculated coupling constants are one order of magnitude smaller than the results from the electronic-structure calculations, regardless of the position of the transition dipole. A similar result was found for HF stretch and HOH bend coupling.101 This finding suggests that the intermolecular coupling between water bend modes is not through transition dipoles, but rather could be, for example, electronic or mechanical. In any case, there is no simple way to include these effects. Meanwhile, we note that there are several sophisticated water models fit from high-level ab-initio calculations, such as the WHBB model developed by Wang et al.102 Compared to accurate electronic-structure calculations such as those at the CCSD(T) level, the computational cost of this model is much cheaper and the accuracy is still excellent. Thus, it is possible to use this model to calculate the bend coupling constants numerically. However, in practice, since we need to calculate numerous coupling terms on the fly, using the WHBB model on a big water cluster at every time step is still expensive. In contrast, to do it on a dimer, on the fly, is very straightforward. How sensitive is the coupling between two molecules in a dimer to the surrounding molecules? To investigate this, we compared the coupling constant between the two bend chromophores of a gas-phase dimer computed by the WHBB

FIG. 3. Comparison between the intermolecular coupling constant k i j [dimer] (in a.u.) for gas-phase dimers obtained by the WHBB model102 and DFT calculations of k i j [cluster] for water clusters. The dashed line is the map result.

model and the value obtained from electronic-structure calculations on the same dimer embedded inside a water cluster. The coupling constant k i j is defined by i j ~ωi j = k i j θ 10 θ 10.

(14)

For the former calculation, k i j for a hydrogen-bonded water dimer that is extracted randomly from the MD simulation is computed by numerical finite differences under the constraints that the center of mass and the HOH bisector of each water molecule are fixed. For the latter calculation, using the same water dimer, we first construct a dimer-centered cluster where the surrounding water molecules are kept fixed. Then, we further keep the center of mass and the HOH bisector of each water molecule in the central dimer unchanged and perform a constrained optimization on this dimer. The coupling constant between these two molecules can be thus obtained by doing

TABLE I. Spectroscopic maps parameterized from TIP4P bulk water simulations. Root-mean-squared deviations (RMSDs) are shown in square brackets. Frequencies and intermolecular couplings are given in units of wavenumber, while the electric field E, coupling constant k i j [dimer], and all other quantities are given in atomic units. HOH bend frequency and angle maps ω 10 = 1581.46 + 2938.51 E ω 21 = 1551.32 + 3147.80 E θ 10 = 0.220 276 − 4.232 17 × 10−5 ω 10 θ 21 = 0.311 43 − 6.001 86 × 10−5 ω 21

[23] [24] [1.0 × 10−4] [1.3 × 10−4]

Dipole and polarizability derivatives and intermolecular coupling map µ ′ = −0.321 63 α ′x x = −0.220 589 α ′y y = 0.544 310 α ′z z = −0.871 378 i θj ω i j = (2.7116 × 105k i j [dimer] + 8.75) θ 10 10

014502-6

Y. Ni and J. L. Skinner

J. Chem. Phys. 143, 014502 (2015)

a frequency calculation on the optimized configuration in Gaussian09. In Fig. 3, we plot the WHBB coupling constant of the dimer versus the ab-initio coupling constant of the cluster. If the points had lain along the diagonal, then we would gather that the coupling constant does not depend on the surrounding molecules. However, the points do not do that. Nonetheless, surprisingly, there is a reasonably good linear relationship (shown by the dashed line) between the WHBB dimer calculation and the cluster calculation. That means that if we know the dimer result, we can obtain the cluster result from this linear map. The full set of spectroscopic maps is summarized in Table I. IV. LINE-SHAPE RESULTS A. Simulation protocols

For the present study, we use a newly developed water model (E3B),103 which explicitly includes three-body interactions in its Hamiltonian, to perform the classical MD simulation. This model has proven to be superior than most of the pair-wise non-polarizable rigid water models in the description of the liquid/vapor interface and the temperature-dependent dynamics of bulk water.25,26,104 MD simulations of liquid H2O within the E3B model are performed in the NVT ensemble with 500 molecules at the experimental densities at 5 ◦C, 25 ◦C, and 60 ◦C using a home-modified GROMACS package version 4.5.5.105 Standard three-dimensional periodic boundary conditions are applied, and electrostatic interactions are calculated using the particle-mesh Ewald summation. For each simulation, the system is maintained at constant temperature by means of the Nosé-Hoover algorithm with the coupling parameter τ = 2 ps. After an equilibration run of 1 ns, a production run lasting 2 ns is performed and the trajectory is sampled every 10 fs for the spectral calculation that follows. The equations of motion are integrated with a 1 fs time step using the SETTLE algorithm for constraints. The simulation of liquid D2O (see below) is applied with the same scheme; only one temperature, 25 ◦C, is sampled. For the simulations of the liquid/vapor interface, the same three temperatures, 5 ◦C, 25 ◦C, and 60 ◦C, are studied in the case of H2O. Using the equilibrated bulk liquid structure obtained above as the initial configuration, the interfacial simulation is performed by tripling the simulation box size in the z direction and otherwise using the same protocol. Similarly, only one temperature, 25 ◦C, is studied in the case of D2O. As discussed in previous studies, the above simulation of the liquid/vapor interface creates a second surface, which also contributes to the calculated spectrum;26,96 in fact, these two surfaces produce cancelling contributions to the SFG spectrum. To overcome this problem, a switching function usually multiplies the transition dipole moments of the chromophores at some time t = 0. This procedure is also used in the present study, and the switching function is given by26  1      f (z) =  (2r c3 + 3r c2 z − z 3)/4r c3   0 

if z > r c if − r c ≤ z ≤ r c , if z < −r c

(15)

where z is the position of the center of mass of a water molecule and z = 0 at the center of the slab. r c is taken to be 4 Å. Therefore, the switching function acts only within 4 Å of the center of the slab, which is well inside the bulk phase. A similar method has been adopted in the earlier study by Nagata et al.,58 although the form of the switching function is different. For the OH stretch, as discussed earlier, one can eliminate the effects of intra- and intermolecular vibrational coupling by studying dilute HOD in D2O. This is obtained experimentally by adding a small amount of H2O to D2O. The equilibrium between H2O, D2O, and HOD readily forms 2 HOD molecules for every H2O molecule. For the bend, there is no intramolecular coupling, but one can eliminate intermolecular coupling by again studying dilute HOD in D2O. In this case, the bend frequency shifts to about 1400 cm−1, making it difficult to compare to the spectroscopy in neat H2O. If instead one increases the amount of H2O added to D2O, at some point one can obtain dilute H2O in a HOD/D2O mixture. This also eliminates the intermolecular vibrational coupling and allows for a direct comparison with neat H2O. Therefore, this is what we pursue here. To calculate the spectra of H2O in HOD/D2O, we assume it is an adequate approximation to consider instead H2O in D2O. In fact, we simulate pure D2O, and for the purposes of the spectroscopy, at each time step in the trajectory, we choose one D2O to be the putative H2O molecule and calculate the electric fields at the molecule from the charges on its surrounding molecules. Then, the transition-frequency and other trajectories can be obtained by using the spectroscopic maps given in Table I. A time average over the trajectories and an ensemble average over each putative H2O molecule result in the final spectral line shape. Because we are simulating the dynamics of D2O, while we are interested in the spectroscopy of H2O, this short-cut (of simulating pure D2O) will only be accurate if the spectroscopy is dominated by static configurations rather than dynamics (motional narrowing), which is the case for the problems of interest (see below). The last parameter needed for the spectral calculation is the lifetime T1 of the HOH bend mode. At room temperature, this value is taken to be 0.17 ps in the case of pure water and 0.38 ps in the case of H2O in D2O, which were determined by experiment.41,106 B. IR line shape for liquid water

The calculated IR line shapes are exhibited in Fig. 4. First, we consider the line shape for neat H2O. Theory (from Eq. (2)) and experiment41 agree well, and both show a peak at 1642 cm−1, which represents a blue-shift of 47 cm−1 from the gas-phase value. This blue-shift comes from hydrogenbonding interactions. In the same figure, we also show experimental results for dilute H2O in HOD/D2O,106 and the calculated results come from Eq. (5) (the intermolecular coupling has been set to zero). Again, agreement is quite good. Theory and experiment show an additional blue-shift of about 10 cm−1, which is due to the absence of intermolecular interactions. This frequency difference is also observed in other experimental measurements.107

014502-7

Y. Ni and J. L. Skinner

J. Chem. Phys. 143, 014502 (2015)

The SFG intensity for ssp polarization (the visible and sum-frequency beams are polarized parallel to the surface and the IR beam is polarized in the scattering plane) is given by40,58 Iss p (ω) ∝| Aeiφ + B χ Rxx z (ω) |2,

(16)

where A is the non-resonant amplitude, φ is the non-resonant phase, and B is the resonant amplitude. This equation can be used to fit the experimental SFG intensity spectrum provided the resonant component χ Rxx z (ω) is known; A, φ, and B are adjustable parameters.58 Alternatively, one can assume a model for χ Rxx z (ω) that involves more adjustable parameters and then fit those as well.40 For our theoretical χ Rxx z (ω), the fitted SFG intensity spectra are shown in the upper panel of Fig. 5. Two fittings are performed based on the different experimental SFG intensity spectral measurements from Nagata et al.58 and Vinaykin and Benderskii.40 Excellent agreement with the experimental result is found in both cases, and both fittings result in a non-resonant phase φ of 40◦. We note that this number is different from the phases determined from other fits (using different models): they obtained −40◦ in Ref. 58 and −57◦ in Ref. 40. In order to compare directly with other theory, and because we would like to interpret these results, in the lower panel of Fig. 5, we show our calculated Im( χ Rxx z (ω)). First considering

FIG. 4. (a) Experimental41,106 and calculated spectra of H2O in D2O and neat H2O at 300 K. All spectra are normalized to have the same peak height. (b) Calculated IR spectra for the HOH bend in neat H2O at 5, 25, 60 ◦C.

We next investigate the spectra for neat H2O at different temperatures and show the results in the lower panel of Fig. 4. By increasing temperature from 5 ◦C to 60 ◦C, the peak positions exhibit a very small (about 1 cm−1) red-shift, which agrees well with the experimental trend.41 This red-shift is caused by the weakening of the hydrogen-bond network between water molecules at high temperature. As the hydrogenbond strength decreases, the restoration force exerted on the HOH bend also decreases, leading to a red-shift of the peak. The line width of the peak also narrows when temperature increases, which is also observed in the experimental IR absorption spectra.41 C. SFG spectrum of the liquid/vapor interface

In recent years, SFG spectra of the HOH bend (in neat H2O) at the liquid/vapor interface have been measured experimentally.40,58 The most informative type of SFG measurement is the heterodyne-detected spectrum, since one can determine the phase, and hence the imaginary part, of the resonant component of the second-order susceptibility of the system, which is analogous to the absorption spectrum. This imaginary part can be calculated directly by using Eq. (7). However, heterodyne-detected SFG spectra of the HOH bend have not yet been obtained experimentally; only the SFG intensity spectra have been measured.40,58

FIG. 5. (a) Experimental40,58 and calculated SFG intensity spectra of the HOH bend at the neat water liquid/vapor interface at 25 ◦C. The theoretical spectra are calculated by fitting to Eq. (16). All spectra are normalized to have the same peak height. (b) Calculated phase-sensitive SFG spectra for the HOH bend at the neat water liquid/vapor interface at 5, 25, 60 ◦C.

014502-8

Y. Ni and J. L. Skinner

J. Chem. Phys. 143, 014502 (2015)

the result at 25 ◦C, we see clearly that there are two features: a negative well at about 1615 cm−1 and a positive peak at about 1670 cm−1. Note that these positions are similar to but not the same as the positions in the theoretical intensity spectrum, indicating the difficulty of trying to interpret the latter. We can compare our result to a classical calculation, using a completely different simulation protocol, by Nagata et al.58 Their result shows similar negative and positive peaks, but at positions of about 1650 and 1730 cm−1, respectively. They were also able to fit the experimental intensity spectrum well, showing the non-uniqueness of fitting with three adjustable parameters (and emphasizing the need to measure the phasesensitive spectrum). In Fig. 5, we also show the temperature dependence of the imaginary part. Both bands show a small red-shift as the temperature increases from 5 ◦C to 60 ◦C, indicating the weakening of the hydrogen-bond network at the surface with increasing thermal energy.

V. MOLECULAR INTERPRETATION OF WATER BEND SPECTROSCOPY

Now that we have demonstrated good agreement between our model and experiment, for both the bulk liquid and the liquid surface, we next use the results of our calculations to provide a molecular interpretation of the spectra. Interpreting the spectra is very difficult in the presence of intermolecular vibrational coupling. Therefore, we consider the isolated H2O chromophore of dilute H2O in D2O. Our molecular interpretation is particularly clean if we focus on the spectral densities instead of the line shapes, because the spectral densities for different hydrogen-bonding classes are strictly additive. The total spectral densities are usually quite similar to the line shapes themselves (due to the relative unimportance of dynamical effects), and so the qualitative interpretation of one informs the other. A. IR spectral density for the bulk liquid

Since the transition dipole for the bend is (approximately) constant, the spectral density for the bend, as stated earlier, is just the distribution of frequencies. This distribution, shown in Fig. 6, is similar to the line shape for the uncoupled H2O in Fig. 4. We next decompose this frequency distribution into contributions from molecules in different hydrogen-bonding environments. As described in the Introduction, we label the environments by the total number of hydrogen bonds involving the molecule of interest and the number of donor hydrogen bonds. We have used a hydrogen-bond definition based on the electronic occupation of the OH σ ∗ orbital on the hydrogenbonding donor molecule, due to partial charge transfer from the lone-pair electrons of the hydrogen-bonding acceptor molecule.108 The main hydrogen-bonding classes in the bulk liquid are 1 N , 2S , 3S , 3 D , 4 D , in increasing order of importance.31 The contributions from each of these classes are shown in Fig. 6. One generally sees that the more hydrogen bonding, the larger the blue-shift. Note also that the distributions for each of the classes overlap significantly.

FIG. 6. (a) Distribution of frequencies for the bend in bulk liquid water at 25 ◦C, with decomposition into sub-distributions for different hydrogenbonding classes. (b) Sub-distributions for the OH-stretch frequencies in bulk liquid water.

It is of interest to compare these frequency distributions to those for OH-stretch spectroscopy. To this end, we theoretically consider dilute HOD or H2O in D2O, focusing on the uncoupled OH stretch. We use our previously developed frequency maps81 to calculate the distribution of OH-stretch frequencies in bulk water, as shown in Fig. 6 panel (b), for the same hydrogen-bonding classes as for the bend. The 3S and 2S distributions each have two features, one for the hydrogenbonded OH and the other for the free OH. The 1 N distribution only has a free-OH peak (since those molecules have no donor hydrogen bonds). These distributions generally show that the more hydrogen bonding, the larger the red-shift, as has been noted many times before. In the Introduction we mentioned that for OH-stretch spectroscopy it was instructive to consider the hydrogenbonding classes of both donor and acceptor molecules. Such an approach might well also be relevant for bend spectroscopy, but here there are possibly two acceptor molecules, and so this kind of analysis gets quite complicated. Therefore, for simplicity, we decided to focus only on the hydrogen-bonding class of the molecule of interest and then do the same for the stretch (so the comparison between stretch and bend is as straightforward as possible).

014502-9

Y. Ni and J. L. Skinner

J. Chem. Phys. 143, 014502 (2015)

The comparison between the bend and stretch frequency distributions suggests that the bend and stretch frequencies are anti-correlated. In fact, this anti-correlation has an interesting history. Some years ago Falk measured the average stretch and bend frequencies for water in a number of different chemical environments. He found that these frequencies were linearly related and satisfied the equation ωb = 1590.6 + 0.2583 (3706 − ω s ),

(17)

where the bend and stretch frequencies, ωb and ω s , respectively, are in cm−1. The average of the symmetric and antisymmetric stretches for gas-phase water is 3706 (and is also very close to 3707, the experimental OH-stretch frequency of gas-phase HOD). The other two numbers were fit parameters. The value 1590.6 is tantalizingly close to the experimental gas-phase bend frequency of 1594.6 (as noted by Falk). This strongly suggests that the equation could be rewritten as ωb − ωbg = m(ω gs − ω s ), ωbg

(18)

ω gs

where and are the gas-phase values. Since m ≈ 0.259 > 0, it means that the average bend and stretch frequencies are anti-correlated. The numerical value indicates that for a given system, a blue-shift (from the gas phase) for the bend is about 1/4 of the red-shift for the stretch. The authors of both of the experimental bend SFG papers40,58 have implicitly invoked this anti-correlation, not for the average frequencies in different chemical environments as above, but rather for the instantaneous frequencies in different local (but chemically identical) environments. A similar idea was considered for the correlation of OH-stretch frequencies with O-O intermolecular distances, comparing different chemical environments with different instantaneous local environments.109 We can test this idea from our theoretical results by forming a scatter plot of bend and stretch frequencies for the same molecule in bulk liquid water. For each molecule, we calculate the bend frequency and average of the two localmode stretch frequencies from our maps and plot this as a point in the bend-stretch frequency plane. For a large and random sampling of molecules, this is shown in Fig. 7. If these points had lain on a line with negative slope, they would be strongly anti-correlated. That is clearly not the case, but they do appear to be weakly anti-correlated. A way to

FIG. 7. Scatter plot of bend versus stretch frequencies, from the maps, for molecules in bulk liquid water. The dashed line is Eq. (18) with m = 0.259.

quantify this iswith the correlation coefficient, defined by r = ⟨δω s δωb ⟩/ ⟨δω2s ⟩⟨δω2b ⟩, where δω refers to the deviation from the average value; in this case, r = −0.252. −1 ≤ r ≤ 1 and so a small negative value is indicative of weak anticorrelation. Returning to the Falk idea, using Eq. (18), and taking the experimental values of ωbg and ω gs , we can fit our data to determine m. The best-fit line is shown on the graph, and the value of m is 0.259, virtually indistinguishable from the value obtained by Falk. Thus, the trend observed by Falk, which on the average the bend blue-shift is about a quarter of the stretch red-shift, is borne out by our results. B. SFG spectral density for the liquid/vapor interface

Moving on to the SFG spectrum for the water surface, in Fig. 8 we show the contributions to the spectral density (Im[ χ Rxx z (ω)] = ⟨a x x (0)m z (0)δ(ω − ω(0))⟩, from Eq. (9)) for the different hydrogen-bonding classes. Note that in this case the contributions can be positive or negative (depending on the sign of m z (0), the projection of the bend transition dipole on the surface normal) and when added together they can cancel, making the intrepretation of the spectrum more difficult. For example, the 3 D spectral density is generally positive, since a typical 3 D molecule has both donor hydrogen bonds pointing

FIG. 8. (a) Main contributions and the total SFG spectral density of the HOH bend for different hydrogen-bond classes at 25 ◦C. (b) Frequency subdistributions for the HOH bend for different hydrogen-bond classes at the water surface.

014502-10

Y. Ni and J. L. Skinner

down to the surface, and hence, the transition dipole (in the direction of the negative of the bisector) is pointing up. Likewise, a 3S molecule with its free-OH pointing straight up has its transition dipole pointing downward, and so its spectral density is generally negative. The other three spectral densities are negative for low frequencies and positive for high frequencies. The total spectral density (sum of all the contributions) is bimodal, with a negative peak at low frequencies and a positive peak at high frequencies. This is similar to the theoretical spectrum in Fig. 5(b), but keep in mind that result includes intermolecular coupling whereas this one does not. From this analysis, we can glean a general understanding of the temperature dependence in Fig. 5(b). As temperature increases, more hydrogen bonds are broken, increasing the populations of the 1 N , 2S , and 3S classes and decreasing the populations of the 3 D and 4 D classes. According to Fig. 8(a), this enhances the negative peak at low frequencies and diminishes the positive peak at high frequencies, in agreement with Fig. 5(b). Nagata et al.58 have proposed that the two peaks come from overlapping but not entirely cancelling negative and positive contributions from those molecules with 0 or 1 donors and those molecules with 2 donors, respectively. Classes 1 N , 2S , and 3S correspond to the former, and classes 3 D and 4 D correspond to the latter, and so we generally agree with this proposal. There are some differences, however, for example, the sum of the 0/1 donor spectral densities (1 N , 2S , and 3S ) leads to a positive contribution between 1650 and 1675 cm−1, contributing to the positive peak (as well as to the negative peak at lower frequencies). For comparison, one can consider the distribution of frequencies for molecules near the surface. Unlike the spectral density, which involves all molecules and uses the lack of centro-symmetry to determine which molecules contribute, here we must make an arbitrary definition of which are the surface molecules. As in previous work,31 we take those that have their oxygen atoms within 2 Å either side of the Gibbs dividing surface (where the density drops to half its bulk value). The different contributions to the frequency distributions for these molecules are shown in Fig. 8 (bottom panel). One sees that for each class of molecules, the distributions are very similar to those of bulk water (see Fig. 6(a)), but of course the relative contributions of each class are quite different. For example, at the surface we see that the 2S molecules make the biggest contribution.31 Finally, we can compare to the theoretical spectral densities for the ssp OH-stretch SFG experiment, which are shown in Fig. 9(a). Note that the total spectral density has three peaks, a sharp positive peak at 3700 cm−1, a negative peak at 3500 cm−1, and a positive peak at 3300 cm−1 (these three peaks match experiments23 on HOD/D2O very well25,26). From looking at the contributions from the different hydrogenbonding classes, one can see here that 4 D and 2S molecules make the biggest contributions. One can see that the sharp high-frequency (free-OH) peak is primarily due to the free OH of the 2S molecules, the negative peak is primarily due to the hydrogen-bonded OH of the 2S molecules, and the low-frequency positive peak is primarily due to the 4 D molecules. (A more sophisticated analysis of this spectrum

J. Chem. Phys. 143, 014502 (2015)

FIG. 9. (a) Main contributions and the total SFG spectral density of the OH stretch for different hydrogen-bond classes. (b) Frequency sub-distributions for the OH stretch for different hydrogen-bond classes at the water surface.

in terms of hydrogen-bond donor-acceptor pairs has recently been presented.31) One can also compare this to the OH-stretch frequency distributions for molecules at the surface (see Fig. 9(b)). The integrated areas of these distributions are of course the same as those for the bend in Fig. 8(b) (since the fractions of each type of molecule at the surface are the same), but now the shapes of the distributions are somewhat different from those for the OH stretch in the bulk, because at the surface, the free-OH bonds are typically pointing out toward the vapor. In any case, one sees that there are sizeable numbers of 3S and 3 D molecules, but unlike for the bend SFG, they do not contribute appreciably to the stretch SFG. This has to do with the different directions of the molecule-fixed transition dipoles: along the OH bond for the stretch, the negative of the bisector for the bend. At the surface, the direction of the OH bond roughly averages to zero for 3S and 3 D molecules, whereas for the bisector, it does not. Therefore, we reach a general conclusion that the bend and stretch SFG spectra are complementary. All the molecular species are at the surface, but the different spectroscopies report on different ones. The bend SFG shows two peaks, with all molecules contributing. The lower frequency peak is primarily from the molecules with fewer hydrogen bonds, and the higher frequency peak is primarily from molecules with more hydrogen bonds. The stretch SFG has three peaks. The free-OH peak in the stretch SFG arises primarily from 2S

014502-11

Y. Ni and J. L. Skinner

molecules. Such molecules do not produce a “free bend,” since they still have one donor. The only analogue of the free-OH for the bend would be from 1 N molecules. Fig. 8(b) shows that these molecules do have the lowest frequencies, and their contribution (see Fig. 8(a)) to the SFG spectrum is negative (as would be expected if both free-OH groups point to the vapor). However, for example, the 2S molecules also have a negative contribution in the same frequency region, and so in this sense, the two higher frequency stretch SFG peaks, with different signs, merge to make one broad bend SFG peak, with the same sign. The lowest frequency stretch SFG peak, primarily from 4 D molecules, in some sense corresponds to the highest frequency bend SFG peak, although in this case it is the 3 D molecules that contribute the most. These arguments are consistent with the anti-correlation between bend and stretch frequencies discussed above.

VI. CONCLUSIONS

In the current study, we extended our mixed quantum/ classical approach used to calculate the IR spectrum of the OH-stretch mode in water to the calculation of the spectrum of the HOH bend mode. Using the ES/MD approach, we first developed ab-initio-based spectroscopic maps for the HOH bend mode. We then calculated the IR line shapes of the H2O/D2O system and neat H2O and found good agreement with the experimental measurements. We also investigated the IR line shapes of the water bend as a function of temperature and successfully reproduced the red-shift of the peak frequency and the narrowing of the line width as temperature increases. We analyzed the IR spectral density in terms of contributions from molecules in different hydrogen-bonding classes. We also calculated SFG spectra of the HOH bend mode at the water liquid/vapor interface. We fit our results to the experimental SFG intensity spectra and found good agreement. We also analyzed the SFG spectral density in terms of contributions from molecules in different hydrogen-bonding classes. The negative peak at lower frequency arises primarily from molecules making one acceptor or one acceptor and one donor hydrogen bonds (leaving two or one free OH, respectively), and the positive peak at higher frequency arises primarily from molecules making one acceptor and one or two donor hydrogen bonds. We also made some effort to compare bend IR and SFG spectroscopies to OH-stretch IR and SFG spectroscopies. We were particularly interested in how the peaks, and their contributions from molecules in different hydrogen-bonding classes, are related in the two sets of spectra. In doing so, we explored the anti-correlation between bend and stretch frequencies, discovered so elegantly 30 yr ago by Falk for water molecules in different chemical environments. We found that the same anti-correlation is weakly evident between bend and stretch frequencies of water molecules in different instantaneous (water) environments. Appreciating this anti-correlation is helpful in understanding the general relationship between bend and stretch IR and SFG spectroscopies. Stretch and bend SFG spectra are found to be complementary, in that the different directions of the molecule-fixed transition-dipole moments

J. Chem. Phys. 143, 014502 (2015)

dictate that the two spectra probe non-identical classes of molecules. Finally, it is probably worth pointing out again that SFG spectra are very difficult to interpret. Part of the problem arises from the contamination of the relevant imaginary part of the susceptibility by the (irrelevant) real part and the nonresonant backgroup in intensity (homodyned) spectra. Part of the problem arises from the complications in any local molecular picture due to intermolecular vibrational coupling, which causes vibrational eigenstates to spread over several or many molecules. These two problems can be overcome with heterodyne detection and isotope dilution (although we do understand that both of these experiments are more difficult). The final problem comes about because SFG spectroscopy involves much cancellation: it is inherently signed, and so different groups of molecules can make contributions with different signs. Moreover, it is a “net” spectroscopy in that even within each group only the net lack of centro-symmetry can be observed. The contributions from different groups of molecules tend to be broad, and their cancellation produces features not assignable solely to a single group of molecules or to single resonances in a phenomenological treatment. With help from molecular simulation, however, SFG spectroscopy remains a magnificently powerful tool for understanding structure and dynamics at surfaces. ACKNOWLEDGMENTS

The authors are grateful for support from NSF Grant No. CHE-1058752. The computational resource for the calculations of spectroscopic maps was supported in part by NSF Grant No. CHE-0840494. 1H.

Bakker and J. L. Skinner, Chem. Rev. 110, 1498 (2010).

2Z. Wang, A. Pakoulev, Y. Pang, and D. D. Dlott, J. Phys. Chem. A 108, 9054

(2004). Palamarev and G. Georgiev, Vib. Spectrosc. 7, 255 (1994). 4D. E. Hare and C. M. Sorensen, J. Chem. Phys. 93, 6954 (1990). 5W. F. Murphy and H. J. Bernstein, J. Phys. Chem. 76, 1147 (1972). 6J. D. Smith, C. D. Cappa, K. R. Wilson, R. C. Cohen, P. L. Geissler, and R. J. Saykally, Proc. Natl. Acad. Sci. U. S. A. 102, 14171 (2005). 7J. B. Asbury, T. Steinel, C. Stromberg, S. A. Corcelli, C. P. Lawrence, J. L. Skinner, and M. D. Fayer, J. Phys. Chem. A 108, 1107 (2004). 8Q. Sun, Chem. Phys. Lett. 568, 90 (2013). 9Q. Wan, L. Spanu, G. A. Galli, and F. Gygi, J. Chem. Theory Comput. 9, 4124 (2013). 10J.-H. Choi and M. Cho, J. Chem. Phys. 138, 174108 (2013). 11P. B. Miranda and Y. R. Shen, J. Phys. Chem. B 103, 3292 (1999). 12M. G. Brown, E. A. Raymond, H. C. Allen, L. F. Scatena, and G. L. Richmond, J. Phys. Chem. A 104, 10220 (2000). 13M. J. Shultz, S. Baldelli, C. Schnitzer, and D. Simonelli, J. Phys. Chem. B 106, 5313 (2002). 14G. L. Richmond, Chem. Rev. 102, 2693 (2002). 15D. Liu, G. Ma, L. M. Levering, and H. C. Allen, J. Phys. Chem. B 108, 2252 (2004). 16Y. R. Shen and V. Ostroverkhov, Chem. Rev. 106, 1140 (2006). 17A. Perry, C. Neipert, B. Space, and P. B. Moore, Chem. Rev. 106, 1234 (2006). 18W. Gan, D. Wu, Z. Zhang, R. Feng, and H. Wang, J. Chem. Phys. 124, 114705 (2006). 19A. Morita and T. Ishiyama, Phys. Chem. Chem. Phys. 10, 5801 (2008). 20X. Fan, X. Chen, L. Yang, P. S. Cremer, and Y. Q. Gao, J. Phys. Chem. B 113, 11672 (2009). 21C. S. Tian, S. J. Byrnes, H.-L. Han, and Y. R. Shen, J. Phys. Chem. Lett. 2, 1946 (2011). 22A. M. Jubb, W. Hua, and H. C. Allen, Acc. Chem. Res. 45, 110 (2011). 3H.

014502-12 23C.-S.

Y. Ni and J. L. Skinner

Tian and Y. R. Shen, J. Am. Chem. Soc. 131, 2790 (2009). S. Tian and Y. R. Shen, Chem. Phys. Lett. 470, 1 (2009). 25P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Am. Chem. Soc. 133, 10360 (2011). 26P. A. Pieniazek, C. J. Tainter, and J. L. Skinner, J. Chem. Phys. 135, 044701 (2011). 27T. Ishiyama and A. Morita, J. Phys. Chem. C 113, 16299 (2009). 28T. Ishiyama and A. Morita, J. Chem. Phys. 131, 244714 (2009). 29S. Nihonyanagi, T. Ishiyama, T.-K. Lee, S. Yamaguchi, M. Bonn, A. Morita, and T. Tahara, J. Am. Chem. Soc. 133, 16875 (2011). 30T. Ishiyama, H. Takahashi, and A. Morita, J. Phys.: Condens. Matter 24, 124107 (2012). 31C. J. Tainter, Y. Ni, L. Shi, and J. L. Skinner, J. Phys. Chem. Lett. 4, 12 (2013). 32B. M. Auer, R. Kumar, J. R. Schmidt, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 104, 14215 (2007). 33M. Falk, Spectrochim. Acta, Part A 40, 43 (1984). 34J. E. Bertie, M. K. Ahmed, and H. H. Eysel, J. Phys. Chem. 93, 2210 (1989). 35J. E. Bertie and Z. Lan, Appl. Spectrosc. 50, 1047 (1996). 36J.-J. Max and C. Chapados, J. Chem. Phys. 131, 184505 (2009). 37D. A. Draegert, N. W. B. Stone, B. Curnutte, and D. Williams, J. Opt. Soc. Am. 56, 64 (1966). 38M. Pavlovic, G. Baranovic, and D. Lovrekovic, Spectrochim. Acta, Part A 47, 897 (1991). 39J. P. Devlin, J. Sadlej, and V. Buch, J. Phys. Chem. A 105, 974 (2001). 40M. Vinaykin and A. V. Benderskii, J. Phys. Chem. Lett. 3, 3348 (2012). 41S. Ashihara, S. Fujioka, and K. Shibuya, Chem. Phys. Lett. 502, 57 (2011). 42S. Ashihara, N. Huse, A. Espagne, E. T. J. Nibbering, and T. Elsaesser, Chem. Phys. Lett. 424, 66 (2006). 43J. Lindner, P. Vöhringer, M. S. Pshenichnikov, D. Cringus, D. A. Wiersma, and M. Mostovoy, Chem. Phys. Lett. 421, 329 (2006). 44P. Bodis, O. F. A. Larsen, and S. Woutersen, J. Phys. Chem. A 109, 5303 (2005). 45N. Huse, S. Ashihara, E. T. J. Nibbering, and T. Elsaesser, Chem. Phys. Lett. 404, 389 (2005). 46O. F. A. Larsen and S. Woutersen, J. Chem. Phys. 121, 12143 (2004). 47F. Ingrosso, R. Rey, T. Elsaesser, and J. T. Hynes, J. Phys. Chem. A 113, 6657 (2009). 48R. Rey and J. T. Hynes, Phys. Chem. Chem. Phys. 14, 6332 (2012). 49A. Bastida, J. Zuniga, A. Requena, and B. Miguel, J. Chem. Phys. 136, 234507 (2012). 50B. Miguel, J. Zuniga, A. Requena, and A. Bastida, J. Phys. Chem. B 118, 9427 (2014). 51J. A. McGuire and Y. R. Shen, Science 313, 1945 (2006). 52M. Smits, A. Ghosh, M. Sterrer, M. Müller, and M. Bonn, Phys. Rev. Lett. 98, 098302 (2007). 53J. R. Scherer, M. K. Go, and S. Kint, J. Phys. Chem. 78, 1304 (1974). 54M. G. Sceats, M. Stavola, and S. A. Rice, J. Chem. Phys. 71, 983 (1979). 55M. Sovago, R. K. Campen, G. W. H. Wurpel, M. Müller, H. J. Bakker, and M. Bonn, Phys. Rev. Lett. 100, 173901 (2008). 56L. D. Marco, K. Ramasesha, and A. Tokmakoff, J. Phys. Chem. B 117, 15319 (2013). 57M. Freda, A. Piluso, A. Santucci, and P. Sassi, Appl. Spectrosc. 59, 1155 (2005). 58Y. Nagata, C.-S. Hsieh, T. Hasegawa, J. Voll, E. H. G. Backus, and M. Bonn, J. Phys. Chem. Lett. 4, 1872 (2013). 59J. Martí, J. A. Padró, and E. Guàrdia, J. Chem. Phys. 105, 639 (1996). 60H. Ahlborn, X. Ji, B. Space, and P. B. Moore, J. Chem. Phys. 111, 10622 (1999). 61A. Perry, C. Neipert, C. R. Kasprzyk, T. Green, B. Space, and P. B. Moore, J. Chem. Phys. 123, 144705 (2005). 62E. Harder, J. D. Eaves, A. Tokmakoff, and B. J. Berne, Proc. Natl. Acad. Sci. U. S. A. 102, 11611 (2005). 63H. Torii, J. Phys. Chem. A 110, 9469 (2006). 64P. K. Mankoo and T. Keyes, J. Chem. Phys. 129, 034504 (2008). 24C.

J. Chem. Phys. 143, 014502 (2015) 65T.

Hasegawa and Y. Tanimura, J. Phys. Chem. B 115, 5545 (2011).

66S. Imoto, S. S. Xantheas, and S. Saito, J. Chem. Phys. 138, 054506 (2013). 67S. Imoto, S. S. Xantheas, and S. Saito, J. Chem. Phys. 139, 044503 (2013). 68L.-P.

Wang, T. L. Head-Gordon, J. W. Ponder, P. Ren, J. D. Chodera, P. K. Eastman, T. J. Martinez, and V. S. Pande, J. Phys. Chem. B 117, 9956 (2013). 69S. Habershon, G. S. Fanourgakis, and D. E. Manolopoulos, J. Chem. Phys. 129, 074501 (2008). 70F. Paesani and G. A. Voth, J. Chem. Phys. 132, 014105 (2010). 71M. Rossi, H. Liu, F. Paesani, J. Bowman, and M. Ceriotti, J. Chem. Phys. 141, 181101 (2014). 72G. R. Medders and F. Paesani, J. Chem. Theory Comput. 11, 1145 (2015). 73P. L. Silvestrelli, M. Bernasconi, and M. Parrinello, Chem. Phys. Lett. 277, 478 (1997). 74W. Chen, M. Sharma, R. Resta, G. Galli, and R. Car, Phys. Rev. B 77, 245114 (2008). 75M. Heyden, J. Sun, S. Funkner, G. Mathias, H. Forbert, M. Havenith, and D. Marx, Proc. Natl. Acad. Sci. U. S. A. 107, 12068 (2010). 76C. Zhang, D. Donadio, and G. Galli, J. Phys. Chem. A 1, 1398 (2010). 77T. Ikeda, J. Chem. Phys. 141, 044501 (2014). 78Y. Nagata, S. Yoshimune, C.-S. Hsieh, J. Hunger, and M. Bonn, Phys. Rev. X 5, 021002 (2015). 79B. M. Auer and J. L. Skinner, J. Chem. Phys. 128, 224511 (2008). 80J. L. Skinner, B. M. Auer, and Y.-S. Lin, Adv. Chem. Phys. 142, 59 (2009). 81S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, J. Chem. Theory Comput. 9, 3109 (2013). 82S. A. Corcelli, C. P. Lawrence, and J. L. Skinner, J. Chem. Phys. 120, 8107 (2004). 83Y.-S. Lin, B. M. Auer, and J. L. Skinner, J. Chem. Phys. 131, 144511 (2009). 84P. A. Pieniazek, Y.-S. Lin, J. Chowdhary, B. M. Ladanyi, and J. L. Skinner, J. Phys. Chem. B 113, 15017 (2009). 85J. L. Skinner, P. A. Pieniazek, and S. M. Gruenbaum, Acc. Chem. Res. 45, 93 (2012). 86C. J. Tainter and J. L. Skinner, J. Chem. Phys. 137, 104304 (2012). 87Y. Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 110, 1992 (2013). 88C. J. Tainter, L. Shi, and J. L. Skinner, J. Chem. Phys. 140, 134503 (2014). 89L. Shi and J. L. Skinner, J. Phys. Chem. B 117, 15536 (2013). 90L. Shi, Y. Ni, S. E. P. Drews, and J. L. Skinner, J. Chem. Phys. 141, 084508 (2014). 91S. Roy, S. M. Gruenbaum, and J. L. Skinner, J. Chem. Phys. 141, 18 (2014). 92S. Roy, S. M. Gruenbaum, and J. L. Skinner, J. Chem. Phys. 141, 22 (2014). 93D. A. McQuarrie, Statistical Mechanics (Harper and Row, New York, 1976). 94F. Li and J. L. Skinner, J. Chem. Phys. 133, 244504 (2010). 95F. Li and J. L. Skinner, J. Chem. Phys. 132, 204505 (2010). 96B. M. Auer and J. L. Skinner, J. Chem. Phys. 129, 214705 (2008). 97B. M. Auer and J. L. Skinner, J. Phys. Chem. B 113, 4125 (2009). 98B. R. Johnson and W. P. Reinhardt, J. Chem. Phys. 85, 4538 (1986). 99D. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1992). 100A. Morita and S. Kato, J. Chem. Phys. 110, 11987 (1999). 101G. Marcotte and P. Ayotte, J. Chem. Phys. 134, 114522 (2011). 102Y. Wang, X. Huang, B. C. Shepler, B. J. Braams, and J. M. Bowman, J. Chem. Phys. 134, 094509 (2011). 103C. J. Tainter, P. A. Pieniazek, Y.-S. Lin, and J. L. Skinner, J. Chem. Phys. 134, 184501 (2011). 104Y. Ni and J. L. Skinner, J. Chem. Phys. 141, 024509 (2014). 105B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. Theory Comput. 4, 435 (2008). 106L. Piatkowski and H. J. Bakker, J. Chem. Phys. 135, 214509 (2011). 107M. Falk, J. Raman Spectrosc. 21, 563 (1990). 108R. Kumar, J. R. Schmidt, and J. L. Skinner, J. Chem. Phys. 126, 204107 (2007). 109C. P. Lawrence and J. L. Skinner, Chem. Phys. Lett. 369, 472 (2003).