THE JOURNAL OF CHEMICAL PHYSICS 142, 212407 (2015)

Molecular dynamics study of two-dimensional sum frequency generation spectra at vapor/water interface Tatsuya Ishiyama,1,a) Akihiro Morita,2,b) and Tahei Tahara3

1

Department of Applied Chemistry, Graduate School of Science and Engineering, University of Toyama, Toyama 930-8555, Japan 2 Department of Chemistry, Graduate School of Science, Tohoku University, Sendai 980-8578, Japan and Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Kyoto 615-8520, Japan 3 Molecular Spectroscopy Laboratory, RIKEN and Ultrafast Spectroscopy Research Team, RIKEN Center for Advanced Photonics (RAP), 2-1 Hirosawa, Wako 351-0198, Japan

(Received 11 January 2015; accepted 18 February 2015; published online 16 March 2015) Two-dimensional heterodyne-detected vibrational sum frequency generation (2D HD-VSFG) spectra at vapor/water interface were studied by molecular dynamics (MD) simulation with a classical flexible and nonpolarizable model. The present model well describes the spectral diffusion of 2D infrared spectrum of bulk water as well as 2D HD-VSFG at the interface. The effect of isotopic dilution on the 2D HD-VSFG was elucidated by comparing the normal (H2O) water and HOD water. We further performed decomposition analysis of 2D HD-VSFG into the hydrogen-bonding and the dangling (or free) OH vibrations, and thereby disentangled the different spectral responses and spectral diffusion in the 2D HD-VSFG. The present MD simulation demonstrated the role of anharmonic coupling between these modes on the cross peak in the 2D HD-VSFG spectrum. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4914299] I. INTRODUCTION

Water surface is of ubiquitous and fundamental importance to many physico-chemical phenomena. Hence, the structure of interfacial water, such as orientation of OH and hydrogen(H)-bonding network at the interface, has been a focus of experimental1,2 and theoretical3,4 studies. Vibrational sum frequency generation (VSFG) spectroscopy5–8 is a powerful experimental tool to probe the interfacial structure and has been applied to analyze the water surface.9,10 From a theoretical viewpoint, we have proposed computational analysis methods of VSFG spectroscopy by molecular dynamics (MD) simulation,11,12 which elaborated the interpretation of the VSFG spectra for aqueous surfaces.13 The computational analysis of the VSFG spectra is subsequently advanced by other groups.14–17 From the experimental side, it is particularly noteworthy that the advent of the phase-sensitive5 or heterodyne-detected (HD)18 VSFG measurement provides the nonlinear susceptibility including its phase, and thereby enables us to distinguish upward and downward OH orientations.19 Close collaboration of the HD VSFG experiment and the MD simulation revealed detailed structural properties of the water surface.20 Apart from the static structural properties, dynamics of water molecules has been also investigated by spectroscopic means and molecular simulation. In the case of bulk liquid water, the dynamics was studied with the time-resolved pumpprobe or photon echo experiments.21 The two-dimensional (2D) infrared (IR) spectroscopy22–26 opened a new avenue to detect the spectral diffusion and chemical exchange dynamics a)Electronic address: [email protected] b)Electronic address: [email protected]

0021-9606/2015/142(21)/212407/13/$30.00

in the femtosecond or picosecond order.27 The MD simulation studies have also aided our understanding about complicated spectral response of bulk water molecules.28–31 The time-resolved techniques of the 2D IR spectroscopy lead to the development of 2D VSFG experiment for investigating the dynamics of the interfacial water. The first time-resolved VSFG study was reported by McGuire and Shen,32 who found great similarities of the H-bonded OH dynamics in the bulk and interfacial water by their pump-probe experiment. Since then, the water dynamics in heterogeneous environments such as air/water, solid/water, and membrane/water interfaces has been investigated by the time-resolved VSFG method.33–37 In recent years, the 2D SFG spectra have been successfully measured by some groups.38–41 The first 2D HD-VSFG spectrum of water surface was reported by Singh et al.42 They observed an off-diagonal cross peak between the dangling and the H-bonding OH frequencies appears promptly after the pump within the time resolution of 0.2 ps. Subsequently Hsieh et al. reported a characteristic slow response of 3500 cm−1 in their 2D HD-VSFG spectrum of water surface,43 though it was not detected in the experiment by Singh et al. Some attempts of MD simulation of the 2D VSFG spectrum were also reported to study the interfacial water dynamics.44,45 In this paper, we extend our computational analysis of VSFG spectroscopy to the time-resolved 2D HD-VSFG spectrum of liquid water using our molecular model. During our computational studies of the VSFG spectroscopy to date, we have developed and refined the molecular models for the SFG calculation.46–48 In particular, the water model for the SFG simulation has been demonstrated to provide accurate and reliable spectra in comparison with the experimental ones.49,50 We confirm the applicability of the water model to describe

142, 212407-1

© 2015 AIP Publishing LLC

212407-2

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

the dynamical properties by calculating the 2D IR spectrum of liquid water. The present computation employed an efficient analytical methodology to calculate the stability matrix,51 the most computationally demanding part, for simulating the 2D HD-VSFG as well as the 2D IR spectra. In the present analysis of the 2D HD-VSFG spectra of water surface, we mainly focus on the mechanism of the above-mentioned cross peak between the frequencies of different OH moieties. For this purpose, we also examined isotopically diluted water, HOD, in addition to the normal water H2O. The present work clearly manifests the importance of the anharmonic coupling besides the dynamic signature of water surface such as the energy transfer. We further investigate the spectral diffusion dynamics of different OH moieties by decomposing the 2D SFG spectrum into the contributions from the dangling and H-bonding OH vibrations. The remainder of this paper is organized as follows. Section II summarizes the methodology of the present computation. The formulation of 2D HD-VSFG spectrum and the approximations introduced here are discussed in

(4) Rabc (τ , τ , τ ) de 3 2 1

Sec. II A, followed by the MD procedures in Sec. II B. The results of the 2D IR spectra are discussed in Sec. III A to confirm the accuracy of calculated dynamical properties of water molecules. Then the 2D HD-VSFG spectra are discussed in Sec. III C. The calculated 2D HD-VSFG spectra are decomposed into two types by introducing a criterion of the H bonds in Sec. III D. Section IV is devoted to conclusion.

II. COMPUTATIONAL METHODOLOGY A. 2D HD-VSFG computation 1. Fourth-order response function

The 2D HD-VSFG spectroscopy is represented with the fourth-order nonlinear response function to the infrared, infrared, infrared and visible lights.44,52 The classical expression of the fourth-order response function is given as44,52

 x−z 1 *    ∂α ab,l (t 3) ∂ µc,k (t 2) + . / = k BT k,l m r (m) f ∂p f ,r (m)(t 2) ∂q f ,r (m)(t 2) ,  x−z     ∂ µ d, j (t 1) ∂ µ˙ e,i (t 0) +   * 1   . / · µ ˙ (t ) µ ˙ (t ) − , d, j 1 e,i 0  ∂q f , s(n)(t 1) ∂p f , s(n)(t 1)   i, j k BT  n f s(n) , - 

(1)

where a, b, c, d, e, f are the components of the Cartesian coordinate (x ∼ z), k B is the Boltzmann constant, and T is the temperature. p f ,r (m) and q f ,r (m) are the f (= x ∼ z) component of the momentum and the Cartesian coordinate, respectively, of rth site of mth molecule. µc,k and α ab,k denote the dipole moment vector and the polarizability tensor of kth molecule, respectively. µ˙ denotes the time derivative of the dipole moment. τ3 = t 3 − t 2, τ2 = t 2 − t 1, and τ1 = t 1 − t 0. Hereafter, we represent the first term of the right hand side of Eq. (1) with RA and the second term with RB, (4),A Rabcde (τ3, τ2, τ1) =

(4),B Rabcde (τ3, τ2, τ1)

  1 *    ∂α ab,l (t 3) ∂ µc,k (t 2) + * 1 . /. µ˙ d, j (t 1) µ˙ e,i (t 0)+/ , k BT k,l m r (m) f ∂p f ,r (m)(t 2) ∂q f ,r (m)(t 2) i, j k BT , -, -

  1 *    ∂α ab,l (t 3) ∂ µc,k (t 2) + *    ∂ µ d, j (t 1) ∂ µ˙ e,i (t 0) + . /. / . =− k BT k,l m r (m) f ∂p f ,r (m)(t 2) ∂q f ,r (m)(t 2) i, j n s(n) f ∂q f , s(n)(t 1) ∂p f , s(n)(t 1) , -, -

2. Molecular properties

The molecular properties of µ, α, and their derivatives in Eq. (1) are represented using the partial charges and charge response kernel (CRK).49 The partial charge at rth site, Q r , is defined by fitting the surrounding electrostatic potential, and the CRK Kr s = ∂Q r /∂Vs is defined by the derivative of the partial charge at the rth site with respect to the electrostatic potential at the sth site.49 These quantities are determined directly by ab initio or density functional theory calculations. Using these quantities, the dipole moment and

(2a)

(2b)

the polarizability of mth molecule are presented as49 µc, m =

site 

Q r (m)qc,r (m),

(3a)

r (m)

α ab, m = −

site site  

Kr s(m)qa,r (m)qb, s(m),

(3b)

r (m) s(m)

where Q r (m) and Kr s(m) are the partial charge at the rth site of mth molecule and the CRK at the sites r and s of the molecule. Equation (3) necessarily reproduces the results of

212407-3

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

the quantum chemical calculations if the partial charges Q are determined by the least square fitting with the constraint of the dipole moment.53 The Q and K are modeled as geometry-dependent quantities as eq

Q r (m) = Q r (m) + eq

mode 

Kr s(m) = Kr s(m) +

∂Q r (m) s t(m), ∂s t(m)

t mode  t

∂Kr s(m) s t(m), ∂s t(m)

···

··· ···

···

Other derivatives can be derived in the same manner.

(4b)

3. Stability matrix

∂p1(t + ∆t) ∂q1(t) .. .

∂p3N (t + ∆t) ∂p3N (t) ∂q1(t + ∆t) ∂p3N (t) .. .

∂p3N (t + ∆t) ∂q1(t) ∂q1(t + ∆t) ∂q1(t) .. .

∂q3N (t + ∆t) ∂p3N (t)

∂q3N (t + ∆t) ∂q1(t)

(7)

where i, k (= 1 ∼ 3N) denote the degrees of freedom and α could be an arbitrary quantity. The stability matrix from t to t + 2∆t can be calculated by the chain rule of the derivatives S(t + 2∆t,t) = S(t + 2∆t,t + ∆t)S(t + ∆t,t)

(8)

or 3N  ∂pi (t + 2∆t)  ∂pi (t + 2∆t) ∂pk (t + ∆t) = ∂p j (t) ∂pk (t + ∆t) ∂p j (t) k=1  ∂pi (t + 2∆t) ∂qk (t + ∆t) + , ∂qk (t + ∆t) ∂p j (t)

(9a)

(5)

Equation (1) contains the derivatives at a same time [∂ µ(t 2)/∂q(t 2) and ∂ µ(t 1)/∂q(t 1)] and those at different times [∂α(t 3)/∂p(t 2) and ∂ µ(t ˙ 0)/∂p(t 1)]. The former derivatives can be analytically calculated using the above molecular model, while the latter requires the calculation of stability matrix S. The stability matrix from t to t + ∆t is defined as

∂p1(t + ∆t) ∂p3N (t) .. .

where 3N is the number of total degrees of freedom. Using the stability matrix S, the derivatives at different times in Eq. (1) are represented, for example,  3N  ∂α(t 3)  ∂α(t 3) ∂pk (t 3) ∂α(t 3) ∂qk (t 3) = + , ∂pi (t 2) k=1 ∂pk (t 3) ∂pi (t 2) ∂qk (t 3) ∂pi (t 2)

  ∂Q s(m) ∂s t(m) ∂ µa, m = Q r (m)δ ab + qa,r (m). ∂qb,r (m) ∂s t(m) ∂qb,r (m) s(m) t

(4a)

where the superscript eq denotes the quantities at equilibrium molecular geometry and s t(m) is the displacement of the tth internal coordinate of mth water molecule.49 K eq, Qeq, ∂K/∂s, ∂Q/∂s are non-empirically determined by ab initio or density functional theory calculations.49 Equations (3) and (4) allow us

∂p1(t + ∆t) *. ∂p1(t) .. .. .. . .. ∂p (t + ∆t) 3N .. ∂p1(t) .. . S(t + ∆t,t) = . ∂q1(t + ∆t) .. ∂p1(t) .. .. .. . .. ∂q (t .. 3N + ∆t) . ∂p1(t) ,

to calculate ∂µ/∂q and ∂α/∂q analytically at an instantaneous configuration. For example,

···

··· ···

···

∂p1(t + ∆t) ∂q3N (t) .. .

+/ // // ∂p3N (t + ∆t) // / ∂q3N (t) // / ∂q1(t + ∆t) // , / ∂q3N (t) // // .. . // ∂q3N (t + ∆t) // / ∂q3N (t) / -

3N  ∂qi (t + 2∆t)  ∂qi (t + 2∆t) ∂pk (t + ∆t) = ∂p j (t) ∂pk (t + ∆t) ∂p j (t) k=1  ∂qi (t + 2∆t) ∂qk (t + ∆t) + , ∂qk (t + ∆t) ∂p j (t) 3N  ∂pi (t + 2∆t)  ∂pi (t + 2∆t) ∂pk (t + ∆t) = ∂q j (t) ∂pk (t + ∆t) ∂q j (t) k=1  ∂pi (t + 2∆t) ∂qk (t + ∆t) , + ∂qk (t + ∆t) ∂q j (t) 3N  ∂qi (t + 2∆t)  ∂qi (t + 2∆t) ∂pk (t + ∆t) = ∂q j (t) ∂pk (t + ∆t) ∂q j (t) k=1  ∂qi (t + 2∆t) ∂qk (t + ∆t) + , ∂qk (t + ∆t) ∂q j (t)

(6)

(9b)

(9c)

(9d)

where i, j, k (= 1 ∼ 3N) denote the degrees of freedom of the whole system. Therefore, the stability matrix for any time interval can be constructed by multiplying the stability matrix for a time step of the MD simulation. We have accordingly implemented the fully analytical calculation method of the

212407-4

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

stability matrix by using the velocity-Verlet algorithm,54 which is much more efficient than the numerical differentiation of different trajectories.30,55 In the velocity-Verlet scheme, the time development of the Cartesian coordinate qi and the momentum pi of ith degree of freedom (i = 1, 2, . . . , 3N) for one time step ∆t are represented as (∆t)2 ∆t pi (t) + Fi (qi (t)), (10a) mi 2mi ∆t ∆t pi (t + ∆t) = pi (t) + Fi (qi (t)) + Fi (qi (t + ∆t)), (10b) 2 2 where mi is the mass for the ith degree of freedom and Fi is the total (intramolecular + intermolecular) force acting on ith coordinate. The elements of the stability matrix in Eq. (6) can be expressed as

distance of r cut = 3.3 Å between the oxygen sites of water was introduced to the calculation of Eq. (12). The distance of 3.3 Å corresponds to the first minimum of the radial distribution function of liquid water. Actually, we confirmed that the calculated results of 2D HD-VSFG are little affected by introducing a cutoff distance larger than 3.3 Å. We note that the calculations of the MD trajectories take account of the long-range forces by the Ewald sum method.54

qi (t + ∆t) = qi (t) +

(∆t)2 ∂Fi (q(t)) qi (t + ∆t) = δi j + q j (t) 2mi ∂q j (t) ∆t Hi j (t), (11a) = δi j + mi pi (t + ∆t) ∆t ∂Fi (q(t)) ∆t ∂Fi (q(t + ∆t)) = + q j (t) 2 ∂q j (t) 2 ∂q j (t) ∆t ∂Fi (q(t)) = 2 ∂q j (t) 3N ∆t  ∂Fi (q(t + ∆t)) ∂qk (t + ∆t) + 2 k=1 ∂qk (t + ∆t) ∂q j (t)   3N  ∆t = Hi j (t) + Hik (t + ∆t) δ k j + Hk j (t) mk k=1 = Hi j (t) + Hi j (t + ∆t) 3N  ∆t + Hik (t + ∆t)Hk j (t), mk k=1 qi (t + ∆t) ∆t = δi j , p j (t) mi

(11b) (11c)

3N ∆t  ∂Fi (q(t + ∆t)) ∂qk (t + ∆t) pi (t + ∆t) = δi j + p j (t) 2 k=1 ∂qk (t + ∆t) ∂p j (t)

(∆t)2 ∂Fi (q(t + ∆t)) 2m j ∂q j (t) ∆t Hi j (t + ∆t), = δi j + mj

= δi j +

(11d)

where Hi j =

∆t ∂Fi (q(t)) . 2 ∂q j (t)

(12)

Hi j is associated to the Hessian, which is the derivative of the force with respect to the Cartesian coordinates. Calculation of H can be carried out analytically, though it is still a computationally demanding task. It would become further demanding by using a polarizable model, as it involves the derivative calculation of the electronic polarization. Therefore, we neglected the effect of electronic polarization in the calculation of Eq. (12). To further reduce the computational cost of Eq. (12), we made use of the short-range character of H in comparison with the force, as the Hessian is a higherorder derivative of energy than the force. Hence, a cutoff

4. Devices for cost reduction

To further reduce the computational cost in the order of magnitude, we introduced the following approximations to the intermolecular contributions in the fourth order response function. These devices allow us to achieve the sufficient statistical sampling with practical computational resources. The effect of these approximations was carefully examined and will be discussed. First, all the intermolecular cross terms in Eq. (1) are ignored, and hence, (4) Rabc (τ , τ , τ ) de 3 2 1  1  *   ∂α ab, m (t 3) ∂ µc, m (t 2) + . / ≈ k BT m r (m) f ∂p f ,r (m)(t 2) ∂q f ,r (m)(t 2) , -

1 · *. µ˙ d, m (t 1) µ˙ e, m (t 0) k BT ,   ∂ µ d, m (t 1) ∂ µ˙ e, m (t 0)  +/ . − ∂q (t ) ∂p (t ) f ,r (m) 1 f ,r (m) 1 r (m) f -

(13)

Second, we were aware that the chain rule calculation in Eq. (9) requires computational cost of O(N 3). The calculation of the stability matrix elements ∂p(t + ∆t)/∂q(t) itself requires the O(N 3) computational cost, as indicated in Eq. (11b). Here, we neglected all the cross terms between different molecules, and thus, only the terms where i, j, k belong to the same molecule are considered in Eqs. (9) and (11b). Third, the effects of intermolecular dielectric coupling on the dipole moment µ and polarizability α were neglected in Eq. (13), so that these molecular properties are represented by those in the isolated environment. This treatment is in accord with the MD simulation using a nonpolarizable model. Consequently, the local field effect56 and the non-Condon effect57 on the spectra are ignored. The above approximations lead to drastic reduction of the computational cost, which is in fact necessary to evaluate the 2D HD-VSFG spectra. However, these approximations may underestimate the effects of the intermolecular coupling in general. On the other hand, the present treatment fully incorporates the effects of the intramolecular coupling and thus should provide a reasonable description for the short-time dynamics. The calculated results presented in this work are corroborated to be consistent to the previous experimental 2D HD-VSFG spectra in the short time, as we discuss in Sec. III. Another basic assumption of the present calculation is the classical treatment of the dynamics, though quantum effects should be significant in the accurate description of the dynamics. However, the quantum effects on the vibrational spectra

212407-5

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

are not apparently substantial in the harmonic oscillators.58 The anharmonic coupling can be also qualitatively described with the classical mechanics, and thus, the present MD simulation will show that the effect of the anharmonic coupling is manifested in the calculated 2D spectra. Quantitative discussion on the spectral intensity and long-time relaxation dynamics should require further accurate treatment based on the quantum mechanics.59 B. MD computation 1. Molecular model

We employed in the present MD simulation a nonpolarizable water model by Marti et al.60 with our modification. The original model of Marti et al. is an extended SPC water to incorporate the intramolecular flexibility and thus has been used in calculating vibrational spectra of water.44,61 However, our preliminary MD simulation showed that the stability matrix calculated with the original model turned out to be fairly unstable within a short time interval for τ1 and τ3. This problem of the original model may be caused by the Morse function for the OH stretching potential, as the second derivative [Eq. (12)] with respect to the OH stretching becomes negative in a long OH distance. Thus, we improved the stability by replacing the Morse function with polynomials for the OH stretching term. Thus, the intramolecular potential function has the following form:49,62 H O

2 uintra =

stability matrix.44 The Ewald sum method is employed for the electrostatic calculations.54 Initially, a slab of water was prepared at the center of the simulation cell, where the interface is normal to the z axis pointing from liquid to vapor. The equilibration run at the temperature 298.15 K was conducted for 30 ps with the velocity scaling method. Then, the 2D SFG computation was carried out in the following scheme in the microcanonical (NVE) ensemble. Figure 1 shows the present scheme to calculate response function Eq. (13), where t 0, t 1, t 2, and t 3 correspond to the time in Eq. (13). The time intervals τ1 and τ3 were sampled up to τmax = 400 MD steps (=400 × 0.305 fs = 122 fs), within which the statistical sampling of the stability matrix is well defined. Regarding the τ2, we examined 4 cases of τ2 values in independent simulations: τ2 = 0.0 ps, 0.1 ps, 0.5 ps, 1.0 ps. The following is the present MD procedure to calculate the response function (see also Fig. 1). Step (a) An MD trajectory is generated for the period of 2τmax + τ2. During the MD trajectory, the molecular Cartesian coordinates q, their momenta p, the polarizabilities α, and the dipole moments µ at each MD time step are saved in an external file. Then, we define the four times during the period illustrated in Fig. 1(a): initial time t 0max, two intermediate times t 1 and t 2, and final time t 3max.

6  [k n (∆r 1)n + k n (∆r 2)n ] + k r r ′∆r 1∆r 2 n=2

+ (kθ /2)(∆r 3)2 + k r θ ∆r 3(∆r 1 + ∆r 2),

(14)

where ∆r 1 and ∆r 2 are the displacements of the two O–H bond lengths from the equilibrium O–H distance and ∆r 3 is that from the equilibrium H–H distance. The equilibrium geometry is set to that of the SPC model.63 Only the first term of r.h.s. is different from the original model by Marti et al., and the potential parameters, k r r ′, kθ , k r θ , are same as Ref. 60. On the other hand, the potential parameters, k n , were reparameterized for the non-polarizable model employed, and we set to k2 = 0.2744, k3 = −0.15, k4 = 0.4, k5 = −9.0, and k6 = 23.0 in the atomic units. The optimization of these parameters were carried out to reproduce the dangling and H-bonding peak positions of the imaginary part of the second order nonlinear susceptibility (see Fig. 5). 2. MD procedure

We enclose 350 water molecules in the simulation cell with dimensions L x × L y × L z = 25 Å × 25 Å × 150 Å. The present system is smaller than the one used for our previous SFG computations (e.g., see Ref. 49) to reduce the computational cost of the 2D spectra. We confirmed that it is still sufficient to describe the interfacial structure and the vibrational spectra of pure water surface.61 The periodic boundary conditions are applied for all three directions. The equation of motions are solved with the velocity-Verlet algorithm54 with a time step of 0.305 fs. Such a small time step (0.305 fs) is necessary for stable calculation of the

FIG. 1. The simulation scheme for the present 2D VSFG calculation. See the text in Sec. II B. (a) MD simulation is carried out from t 0 to t 3. (b) TCF is calculated in forward and backward directions. (c) MD simulation is carried out from t 3 to t 3 + ∆t, and calculate a new TCF.

212407-6

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

These times are defined so that t 1 − t 0max = τmax, t 2 − t 1 = τ2, and t 3max − t 2 = τmax are satisfied. Step (b) After the MD simulation at t 3max, the stability matrices from t 1 to t 0 (backward) and from t 2 to t 3 (forward) are calculated up to the time interval τmax using the stored trajectory data. Then, we calculate the derivative quantities included in Eq. (13) and thereby construct the statistical samples of the time TCF in Eq. (13). The above procedure of step (b) is associated to the MD simulation step at t 3max. We call this set of procedure one cycle. Note that t 0 (t 3) varies from t 1 to t 0max (from t 2 to t 3max), and accordingly, one cycle yields the samples of R(4)(τ3, τ2, τ1) with 0 ≤ τ1(τ3) ≤ τmax for a fixed value of τ2. Step (c) To proceed to the second cycle, we move forward the MD simulation from t 3max by the time step ∆t and save the data in the file. Then, all the times of t 0max, t 1, t 2, and t 3max are shifted by ∆t, and the procedure of step (b) is performed using the new times. This cycle is repeated until the statistical average of the response function converges. The converged response functions are Fourier-Laplacetransformed to get the rephasing and nonrephasing signals45 as  ∞ RE R (ω3, τ2,ω1) = Re dτ3eiω3τ3 0



∞

×

 dτ1R(τ3, τ2, τ1)e−iω1τ1 , (15a)

0 ∞

 RNR(ω3, τ2,ω1) = Re

dτ3eiω3τ3

0

 ×

∞

 dτ1R(τ3, τ2, τ1)eiω1τ1 .

(15b)

0

In the numerical calculation of Eq. (15), the Fourier-Laplace transformation was performed up to τmax. The total 2D SFG spectrum was obtained as the sum of them, RT = RRE + RNR. The statistical sampling of the 2D spectrum at each τ2 was taken using 32 independent MD trajectories, which were generated with parallel computers. Each trajectory was utilized to sample 9000 cycles, and thus, the total statistical samples amount to 9000 × 32 = 288 000 cycles for the 2D spectrum at each τ2. All the response functions calculated in the present work are the RTx x z z z component.

C. 2D IR spectrum

We also calculated the third order response function R(3) and the 2D IR spectrum to assess the bulk dynamics of the present MD system.64 Most of the calculation methodology are similar to the above-mentioned procedure for the 2D SFG calculation. In the case of the 2D IR spectrum, the response function was calculated by the following formula instead of Eq. (13):

R(3) x x x x (τ3, τ2, τ1)  x−z ∂ µ x (t 2) + 1  *   ∂ µ x (t 3) . / ≈ k BT m r (m) f ∂p f ,r (m)(t 2) ∂q f ,r (m)(t 2) , 1 · *. µ˙ x (t 1) µ˙ x (t 0) k BT ,  x−z  ∂ µ x (t 1) ∂ µ˙ x (t 0) + / . − (16) ∂q f ,r (m)(t 1) ∂p f ,r (m)(t 1) r (m) f Here, we also denote the first term of r.h.s. of Eq. (16) as RA, the second term as RB, in the same manner as in Eq. (2). To calculate Eq. (16), we prepared a bulk system consisting of 350 water molecules contained in a simulation box with the dimensions of L x × L y × L z = 21.8 Å × 21.8 Å × 21.8 Å to reproduce the bulk water density. The other MD conditions were the same as in the 2D SFG calculation described in Sec. II B. The rephasing and nonrephasing spectra of the 2D IR were calculated by the Fourier-Laplace transformation of R(3) in Eq. (15), and the total 2D IR spectrum was given with the sum of them. III. RESULTS AND DISCUSSION A. 2D IR spectra

First, we discuss the calculated 2D IR spectra for bulk water to examine the accuracy of the present simulation. Figure 2 shows the various terms of the calculated 2D IR spectra of the bulk H2O at the waiting time τ2 = 0 ps. Panel (a) is the R(3),A term, panel (b) is the R(3),B term, panel (c) is the rephasing signal R(3),RE, panel (d) is the nonrephasing signal R(3),NR, and panel (e) is the total 2D IR spectrum R(3),T. Note that R(3),T = R(3),A + R(3),B = R(3),RE + R(3),NR holds for the spectra. Panels (a) and (b) show that RA and RB have very similar shapes with opposite sign. Thus, significant cancellation of RA and RB occurs to yield the total 2D IR spectrum R(3),T in panel (e), which is consistent to the previous work.55 The remainder of the cancellation is attributed to the anharmonicity of the system.65 In fact, we confirmed by our MD simulation and discussed in the supplementary material66 that for a completely harmonic system, the RA and RB terms fully cancel each other, and thus, no 2D IR spectrum is observed. These statements are also valid for the 2D SFG spectra discussed later. In Figs. 2(c) and 2(d), R(3),RE is diagonally elongated, and R(3),NR is anti-diagonally elongated.67 These features come from the phase-twisted mechanism of the rephasing and nonrephasing spectra.68 The sum of the rephasing and nonrephasing signals cancels the dispersive components, yielding a purely absorptive 2D IR spectrum shown in Fig. 2(e). The bleach band (v = 0 → 1 transition) of the 2D IR spectrum appears as a blue-colored region in Fig. 2(e). The red-colored band in Fig. 2(e) is assigned to an anharmonically shifted v = 1 → 2 absorption, so-called hot band. The experimentally observed R(3),RE, R(3),NR, and R(3),T spectra67 are shown in Fig. 3. We find that the calculated spectra well reproduce the experimental features of the 2D

212407-7

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

FIG. 3. The experimental 2D IR spectra of 1% solution of HOD in D2O water system67 for the rephasing, the nonrephasing, and the total signals at τ 2 = 0.1 ps. (Reprinted with permission from J. Chem. Phys. 125, 194521 (2006). Copyright 2006 American Institute of Physics.) (Note that the coloring in the experimental spectra is different from that in the present calculations in Fig. 2 or 4.)

FIG. 2. The calculated 2D IR spectra of the bulk H2O water system for (a) R A term, (b) R B term, (c) the rephasing, (d) the nonrephasing, and (e) the total signals at the waiting time τ 2 = 0 ps.

IR spectra, including the bleach and hot bands and the skewed feature. The time evolution of the calculated total 2D IR spectra R(3),T with the waiting period τ2 is shown in Fig. 4(a) and compared with the experimental one in Fig. 4(b). The partial

2D IR spectra R(3),A are also shown. As τ2 evolves, the diagonal elongation in both R(3),T and R(3),A decays and eventually disappears, which is indicative of the spectral diffusion. This result suggests that the spectral response of the partial response function R(3),A (and R(3),B) reflects the essential features of the total response function R(3),T. The experimental spectra in panel (b) exhibit quite analogous evolution with a comparable time scale to the present calculation. We note that the experimental spectra were measured for the isotopically diluted HOD, where the interand intramolecular vibrational couplings are eliminated. In fact we also examined the 2D IR of the HOD water by the MD simulation and found that the decay of the elongation for HOD water is quite similar with H2O. Qualitative agreement of the experimental 2D IR spectra to the present calculation implies that the time evolution of the spectral diffusion well reflects the dynamics of liquid water. B. 1D VSFG spectra

Before discussing the 2D VSFG spectrum, we examine the steady-state HD-VSFG spectrum (1D HD-VSFG spectrum) of water surface calculated by the present nonpolarizable water model. In our previous studies of 1D VSFG spectra, we employed the polarizable water models

212407-8

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

FIG. 4. (a) The calculated partial and (3),T (3),A total 2D IR spectra, R x x x x and R x x x x of the bulk H2O water system. (b) The experimental 2D IR spectra of 1% solution of HOD in D2O water system.67 (Reprinted with permission from J. Chem. Phys. 125, 194521 (2006). Copyright 2006 American Institute of Physics.)

for the VSFG calculation49,69 with fully taking account of the intermolecular polarization couplings and demonstrated that these calculations are quite accurate to reproduce the 1D HD-VSFG spectra. However, we employed a non-polarizable model to reduce the computational cost of 2D HD-VSFG calculation as we explained in Sec. II B. Therefore, we need to examine the 1D spectrum by using the present water model. The 1D HD-VSFG was calculated by the time correlation function (TCF) formalism56,70 with the present model, where the instantaneous polarizability and the dipole moment of the constituent molecules were calculated by

Eqs. (3) and (4), and all the intermolecular cross terms of the TCF, α i (t)µ j (0) (i , j), were ignored so as to be consistent to the calculation of 2D SFG. The calculations were performed for pure H2O and pure HOD water. Figure 5 shows the calculated imaginary part of the second-order nonlinear susceptibility Im[ χ(2) x x z ]. This xxx component corresponds to the SSP polarization combination of SFG measurement, which consists of S-polarized SFG, S-polarized visible, and P-polarized infrared lights. In the OH stretching vibrations for 3100-3800 cm−1, one can see the positive band of the dangling OH at about 3700 cm−1 and

212407-9

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

FIG. 5. The calculated imaginary part of the second-order nonlinear suscep(2) tibility χ x x z of the present non-polarizable model for full H2O (blue) and full HOD (red) water. The inset is the calculated imaginary part of H2O by the polarizable model.20

the negative band of the H-bonding OH below 3600 cm−1. These signs of the imaginary χ(2) reflect the upward and downward directions of OH, respectively.13 These qualitative features of the signs are commonly seen in H2O water and HOD water, though the amplitude of H2O is almost doubled due to the difference in the OH number density. The H2O water shows more pronounced tail of Im[ χ(2) x x z ] in the low frequency region. We note that the present MD calculation can include the coupled vibrational dynamics of liquid H2O and hence reproduces the low-frequency vibrations in the MD trajectories, while the intermolecular coupling of polarization at each time step is neglected. Comparing the present result of Im[ χ(2) x x z ] with that of the previous polarizable simulation (inset of Fig. 5),20 one can see following two discrepancies. First, the positive feature at the low frequency region below 3150 cm−1 is missing in the present simulation. This is fully in accord with the present simulation that neglects the intermolecular coupling of polarization in the liquid environment, since the coupling of polarization has been proved to be critical for the low-frequency positive feature.50,69 Second, the calculated negative amplitude of the H-bonding band is relatively weak compared to the positive amplitude of the dangling peak. The underestimated amplitude comes from the fact that the present nonpolarizable model neglects the enhanced transition dipole of OH stretching by the H-bonds, which is relevant to the non-Condon effect.57 In the following discussion on the 2D SFG, the present model is applicable to discuss qualitative features of the 2D VSFG spectra, including the decomposition of the dangling and Hbonding OH modes, though the low-frequency tail and the quantitative amplitude in the H-bonding region may be beyond the current scope. C. 2D HD-VSFG spectra

Next, we discuss the 2D HD-VSFG spectra of water surface of pure H2O and pure HOD molecules. An observation of difference between two cases enables us to discuss a role of vibrational couplings. The calculated 2D HD-VSFG spectra RTx x z z z of the H2O and HOD water surfaces are shown in Figs. 6 and 7, respectively. First, we focus on the spectra at the waiting time τ2 = 0.0 ps (the top panels of Figs. 6 and 7). In either panel, the strong negative blue-colored band is seen in the upper

(4),T

FIG. 6. The calculated 2D HD-VSFG spectra R x x z z z for H2O case.

right region, and a positive red-colored band is present below the blue band. The former is assigned to the bleach (0-1 transition) band of the dangling OH mode, and the latter is a mixture of the hot band of the dangling OH and the bleach band of the H-bonding OH. Note that the bleach bands along the diagonal line of ν˜1 ≈ ν˜3 should change their signs between the dangling and H-bonding frequency regions, as indicated by the signs of Fig. 5. A noticeable difference between H2O and HOD is observed in the upper left region; the negative band extends to the upper left region in the panel of H2O, while such feature is not seen in the panel of HOD. This feature of H2O is consistent to the experimentally observed 2D SFG spectrum of H2O.42 This means that the pumping of OH vibration at the H-bonding region ( ν˜1 = 3200–3400 cm−1) immediately (τ2 = 0) decreases the signal of the dangling OH stretching

212407-10

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

The cross peak at H2O surface gradually disappears as τ2 evolves, and the spectra of H2O and HOD surfaces resemble each other at about τ2 = 1.0 ps due to loss of memory of vibration.23 Compared to the experimental 2D SFG spectrum,42 the calculated spectra obviously tend to emphasize the dangling OH response in comparison with the H-bonding response. This tendency is attributed to the missing polarization effect in the present model, as discussed in Fig. 5. To obtain further insight into the dangling and H-bonding OH responses, we carried out decomposition analysis in the following. D. Decomposition of 2D HD-VSFG spectra for HOD water

We have argued above that the hot band of the dangling OH and the bleach band of the H-bonding OH are superposed in the total 2D HD-VSFG spectrum, which complicates the assignment and interpretation of 2D HD-VSFG spectrum. Analysis of OH vibrational spectra of water is generally complicated by various overlap and couplings, and the HD isotope dilution is often invoked to simplify the spectral assignment. However, the H–D dilution is not sufficient to remove the spectral overlap of the dangling and H-bonding signals in the 2D HD-VSFG spectra as discussed above. Therefore, we decompose the 2D HD-VSFG spectra of HOD into the dangling and H-bonding signals by the following scheme. 1. H-bonding OH and dangling OH

(4),T

FIG. 7. The calculated 2D HD-VSFG spectra R x x z z z for HOD case.

mode (ν˜3 ≈ 3700 cm−1) at the H2O surface. This cross peak could be interpreted with an immediate intramolecular energy transfer from the H-bonding to the dangling OH vibrations and/or the anharmonic coupling between these modes.42 Although the opposite energy transfers from the dangling ( ν˜1 ≈ 3700 cm−1) to the H-bonding ( ν˜3 = 3200–3400 cm−1), OH should make another cross peak in the lower right region, such feature is obscured by the overlap of the hot band of the dangling OH. Another difference is seen in the lower left corner; the 2D SFG spectrum of HOD shows a blue region in contrast to that of H2O. This blue band of HOD is attributed to the hot band of the H-bonding OH. This feature is not clear in the H2O spectrum, since the H2O water exhibits more pronounced low-frequency tail of the H-bonding OH fundamental vibrations, as shown in Fig. 5.

The HOD molecules at the surface can be classified into two components: those having hydrogen bonding OH and having dangling OH. In the present MD simulation, the decomposition of the 2D HD-VSFG spectra into two components is straightforward, because the response function Eq. (13) consists of the sum of the individual molecules. To distinguish the dangling and H-bonding OH, the following geometric criteria for the H-bond are employed;49 the H · · · O distance is less than 2.5 Å and the ∠O–H · · · O angle is greater than 140◦. Such decomposition of 2D VSFG might become invalid if the conformational change took place during the waiting time τ2. In fact, we examined the conformation about the H bond with the above criteria during the waiting time τ2 and confirmed that the conformational exchange can be neglected within 1 ps. This result is in accord with the previous study of the H-bond dynamics at water surface45 and validates the above decomposition. Figures 8 and 9 show the decomposed 2D HD-VSFG spectra of the HOD into the H-bonding and dangling OH, respectively. The sum of the two spectra at τ2 = 0 ps (top panels in Figs. 8 and 9) approximately corresponds to the full spectrum of HOD (top panel in Fig. 7). Comparing the two Figures of 8 and 9, we call for attention the different scale of amplitude. This difference is pertinent to the tendency that the non-polarizable simulation relatively underestimates the Hbonding component, as discussed in Sec. III B. In the actual experiment or polarizable simulation, the H-bonding feature should be enhanced because of the non-Condon effect.57 The

212407-11

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

(4),T

FIG. 8. The calculated 2D HD-VSFG spectra R x x z z z only for HOD water with a bonding hydrogen (bonding H).

present decomposition analysis rather allows for inspecting the two components separately in the present MD results. In Fig. 8, the H-bonding OH component at τ2 = 0 ps clearly shows that the red-colored bleach band is elongated along the diagonal line and is accompanied with the bluecolored hot band below. This elongation is indicative of the inhomogeneous broadening, similar to the bulk 2D IR spectrum67 in Fig. 4. This feature of diagonal elongation gradually decays on the time scale of a few hundred femtoseconds. This decaying behavior and the time scale are very similar to those of the 2D IR spectra in Fig. 4, indicating that the spectral diffusion dynamics of H-bonding OH at the water surface is similar to that in the bulk water which is consistent to the pump-probe experiment by McGuire and Shen.32 On

(4),T

FIG. 9. The calculated 2D HD-VSFG spectra R x x z z z only for HOD water with a non-bonding hydrogen (dangling H).

the other hand, the dangling OH component in Fig. 9 exhibits no feature of the diagonal elongation in contrast to the Hbonding component. This is a reasonable result because the dangling OH is located in a more isolated, homogeneous environment than the H-bonding OH. 2. Anharmonic coupling of OH and OD

Next, we discuss the effect of intramolecular anharmonic coupling of OH and OD in the HOD molecule. In the case of dangling OH component discussed in Fig. 9, the HOD molecules at the surface also have OD bonds, which are mostly hydrogen-bonded. To investigate the vibrational

212407-12

Ishiyama, Morita, and Tahara

J. Chem. Phys. 142, 212407 (2015)

by eliminating all the anharmonic coupling, as discussed in the supplementary material.66 These cross peaks are not a consequence of dynamics of molecules, including energy transfer between different vibrational modes, but due to an intrinsic nature of the 2D spectra as long as the potential of the observed system contains anharmonicity.

IV. SUMMARY

(4),T

FIG. 10. (a) The calculated 2D HD-VSFG spectra R x x z z z only for HOD water with a non-bonding hydrogen (dangling H) at τ 2 = 0 ps, which is the same as the top panel of Fig. 9. (b) Energy level scheme of HOD water with a non-bonding hydrogen (dangling H), where the colors of the arrows are in accord with those of the bands in panel (a).

coupling in this case, the 2D HD-VSFG spectrum of Fig. 9 is displayed in a wide frequency range which encompasses both the OH and OD stretching. Figure 10(a) shows the extended spectrum in the case of dangling OH at τ2 = 0 ps from 2000 to 4000 cm−1 region, where the spectrum in the upper right (3100-3900 cm−1) region is identical to the top panel of Fig. 9. The 2D VSFG spectrum in Figure 10(a) clearly manifests the bleach bands, the hot bands, and cross bands of the dangling OH and hydrogen-bonded OD vibrational modes. All the bands can be assigned with the level scheme of two coupled oscillators in Figure 10(b),21,65 where the two numbers in the ket states (e.g., |00⟩, |10⟩, . . . ) indicates the quantum numbers of the OH and OD vibrational modes, respectively. The transitions ①-⑧ shown in the level scheme (b) manifest themselves in the 2D HD-VSFG spectrum in panel (a). The transitions ⑧ and ④ (② and ⑥) correspond to the bleach and the stimulated emission of OH (OD) vibrational mode, respectively, and appear in panel (a) along the diagonal line. We notice that the bands of the OH and OD vibrations have opposite signs, as panel (a) shows the ⑧+④ band in blue the ②+⑥ band in red. This sign relation of the OH and OD bands for the 2D HD-VSFG is different from that of the conventional 2D IR spectra of coupled oscillators,21,65 because the dangling OH and hydrogen-bonded OD pose opposite directions with respect to the interface. The transitions ③ and ⑤ are assigned to the excited state absorption (hot band) of OH and OD. The cross peaks in Fig. 10(a) are assigned to the transitions ①, ②, ⑦, and ⑧, due to the anharmonic coupling of OH and OD modes. We note that without the anharmonic coupling, the band splitting of ① and ② and that of ⑦ and ⑧ should not occur, and thus, these cross peaks would cancel each other. We also confirmed the disappearance of the cross peaks

In this work, we extended our computational method of the vibrational SFG simulation to the two-dimensional VSFG of water interface. Evaluation of the fourth-order response function was facilitated by using the efficient analytical computation of the stability matrix. We employed a nonpolarizable model and neglected the intermolecular couplings to reduce the computational cost, though the calculated 2D IR spectra for the bulk water and the 2D VSFG spectra for the water surface well reproduced the qualitative features of band shapes and relaxation of the corresponding experimental spectra. The complicated 2D VSFG spectra were analyzed by decomposing them into the H-bonded and dangling OH components, which greatly clarified the coupling mechanism of the 2D spectra and the relaxation dynamics of each component. In the isotopically diluted water, the present MD simulation predicted the cross bands of OH and OD modes, which clearly manifest the anharmonic coupling of these modes. ACKNOWLEDGMENTS

We are grateful to Professor Shoichi Yamaguchi, Dr. Satoshi Nihonyanagi, and Dr. Ken-ichi Inoue for stimulating discussions. This work was supported by the Grantin-Aid for Scientific Research on Innovative Areas “Soft Molecular Systems,” Grant-in-Aid for Scientific Research (A) (Grant No. 24245006), and the Next Generation Supercomputer Project of MEXT, Japan. 1P.

B. Petersen and R. J. Saykally, Annu. Rev. Phys. Chem. 57, 333 (2006). Davidovits, C. E. Kolb, L. R. Williams, J. T. Jayne, and D. R. Worsnop, Nature 491, 582 (2012). 3D. Chandler, Nature 437, 640 (2005). 4P. Jungwirth and D. J. Tobias, Chem. Rev. 106, 1259 (2006). 5C. S. Tian and Y. R. Shen, Surf. Sci. Rep. 69, 105 (2014). 6G. L. Richmond, Chem. Rev. 102, 2693 (2002). 7S. Gopalakrishnan, D. Liu, H. C. Allen, M. Kuo, and M. J. Shultz, Chem. Rev. 106, 1155 (2006). 8P. L. Geissler, Annu. Rev. Phys. Chem. 64, 317 (2013). 9Y. R. Shen and V. Ostroverkhov, Chem. Rev. 106, 1140 (2006). 10M. Sovago, R. K. Campen, H. J. Bakker, and M. Bonn, Chem. Phys. Lett. 470, 7 (2009). 11A. Morita and J. T. Hynes, Chem. Phys. 258, 371 (2000). 12A. Morita and J. T. Hynes, J. Phys. Chem. B 106, 673 (2002). 13T. Ishiyama, T. Imamura, and A. Morita, Chem. Rev. 114, 8447 (2014). 14J. L. Skinner, P. A. Pieniazek, and S. M. Gruenbaum, Acc. Chem. Res. 45, 93 (2012). 15Y. Nagata, R. Pool, E. H. G. Backus, and M. Bonn, Phys. Rev. Lett. 109, 226101 (2012). 16M. Sulpizi, M. Salanne, M. Sprik, and M. Gaigeot, J. Phys. Chem. Lett. 4, 83 (2013). 17S. Roy and D. Hore, J. Phys. Chem. C 116, 22867 (2012). 18S. Nihonyanagi, J. A. Mondal, S. Yamaguchi, and T. Tahara, Annu. Rev. Phys. Chem. 64, 579 (2013). 19S. Nihonyanagi, S. Yamaguchi, and T. Tahara, J. Chem. Phys. 130, 204704 (2009). 2P.

212407-13 20S.

Ishiyama, Morita, and Tahara

Nihonyanagi, T. Ishiyama, T. Lee, S. Yamaguchi, M. Bonn, A. Morita, and T. Tahara, J. Am. Chem. Soc. 133, 16875 (2011). 21V. Cervetto, J. Helbing, J. Bredenbeck, and P. Hamm, J. Chem. Phys. 121, 5935 (2004). 22M. D. Fayer, Ultrafast Infrared Vibrational Spectroscopy (CRC Press, Boca Raton, 2013). 23S. T. Roberts, K. Ramasesha, and A. Tokmakoff, Acc. Chem. Res. 42, 1239 (2009). 24M. L. Cowan, B. D. Bruner, N. Huse, J. R. Dwyer, B. Chugh, E. T. J. Nibbering, T. Elsaesser, and R. J. D. Miller, Nature 434, 199 (2005). 25T. E. Erik and T. J. Nibbering, Chem. Rev. 104, 1887 (2004). 26Y. S. Kim and R. M. Hochstrasser, Proc. Natl. Acad. Sci. U. S. A. 102, 11185 (2008). 27M. D. Fayer, Watching Ultrafast Molecular Motions with 2D IR Chemical Exchange Spectroscopy (World Scientific Publishing, Singapore, 2011). 28H. J. Bakker and J. L. Skinner, Chem. Rev. 110, 1498 (2010). 29T. Yagasaki and S. Saito, Annu. Rev. Phys. Chem. 64, 55 (2013). 30T. Hasegawa and Y. Tanimura, J. Chem. Phys. 128, 064511 (2008). 31G. Stirnemann, P. J. Rossky, J. T. Hynes, and D. Laage, Faraday Discuss. 146, 263 (2010). 32J. A. McGuire and Y. R. Shen, Science 313, 1945 (2006). 33A. Ghosh, M. Smits, J. Bredenbeck, and M. Bonn, J. Am. Chem. Soc. 129, 9608 (2007). 34A. Eftekhari-Bafrooei and E. Borguet, J. Am. Chem. Soc. 131, 12034 (2009). 35S. Roy, S. M. Gruenbaum, and J. L. Skinner, J. Chem. Phys. 141, 22D505 (2014). 36C.-S. Hsieh, R. K. Campen, A. C. V. Verde, P. Bolhuis, H.-K. Nienhuys, and M. Bonn, Phys. Rev. Lett. 107, 116102 (2011). 37C.-S. Hsieh, R. K. Campen, M. Okuno, E. H. G. Backus, Y. Nagata, and M. Bonn, Proc. Natl. Acad. Sci. U. S. A. 110, 18780 (2013). 38W. Xiong, J. E. Laaser, R. D. Mehlenbacher, and M. T. Zanni, Proc. Natl. Acad. Sci. U. S. A. 108, 20902 (2011). 39Z. Zhang, L. Piatkowski, H. J. Bakker, and M. Bonn, Nat. Chem. 3, 888 (2011). 40P. Singh, S. Nihonyanagi, S. Yamaguchi, and T. Tahara, J. Chem. Phys. 137, 094706 (2012). 41Z. Zhang, L. Piatkowski, H. J. Bakker, and M. Bonn, J. Chem. Phys. 135, 021101 (2011). 42P. Singh, S. Nihonyanagi, S. Yamaguchi, and T. Tahara, J. Chem. Phys. 139, 161101 (2013). 43C.-S. Hsieh, M. Okuno, J. Hunger, E. H. G. Backus, Y. Nagata, and M. Bonn, Angew. Chem., Int. Ed. 53, 8146 (2014).

J. Chem. Phys. 142, 212407 (2015) 44Y.

Nagata and S. Mukamel, J. Am. Chem. Soc. 133, 3276 (2011). Ni, S. M. Gruenbaum, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 110, 1992 (2013). 46T. Ishiyama, V. V. Sokolov, and A. Morita, J. Chem. Phys. 134, 024509 (2011). 47T. Ishiyama and A. Morita, J. Phys. Chem. C 115, 13704 (2011). 48T. Kawaguchi, K. Shiratori, Y. Henmi, T. Ishiyama, and A. Morita, J. Phys. Chem. C 116, 13169 (2012). 49T. Ishiyama and A. Morita, J. Chem. Phys. 131, 244714 (2009). 50T. Ishiyama, H. Takahashi, and A. Morita, Phys. Rev. B 86, 035408 (2012). 51S. Sakaguchi, T. Ishiyama, and A. Morita, J. Chem. Phys. 140, 144109 (2014). 52S. Mukamel, Principles of Nonlinear Optical Spectroscopy (Oxford University Press, Oxford, 1995). 53T. Ishida and A. Morita, J. Chem. Phys. 125, 074112 (2006). 54M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987). 55J. Jeon and M. Cho, New J. Phys. 12, 065001 (2010). 56A. Morita and T. Ishiyama, Phys. Chem. Chem. Phys. 10, 5801 (2008). 57J. R. Schmidt, S. A. Corcelli, and J. L. Skinner, J. Chem. Phys. 123, 044513 (2005). 58J. S. Bader and B. J. Berne, J. Chem. Phys. 100, 8359 (1994). 59G. Stock, Phys. Rev. Lett. 102, 118301 (2009). 60J. Marti, J. A. Padro, and E. Guardia, J. Mol. Liq. 62, 17 (1994). 61A. Morita, J. Phys. Chem. B 110, 3158 (2006). 62T. Ishiyama and A. Morita, J. Phys. Chem. C 111, 721 (2007). 63H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, “Interaction models for water in relation to protein hydration,” in Intermolecular Forces, edited by B. Pullman (Reidel, Dordrecht, 1981). 64J. R. Schmidt, S. T. Roberts, J. J. Loparo, A. Tokmakoff, M. D. Fayer, and J. Skinner, Chem. Phys. 341, 143 (2007). 65P. Hamm and M. T. Zanni, Concepts and Methods of 2D Infrared Spectroscopy (Cambridge University Press, Cambridge, 2011). 66See supplementary material at http://dx.doi.org/10.1063/1.4914299 for decomposed spectra of 2D response functions and the case of a harmonic system. 67J. Loparo, S. Roberts, and A. Tokmakoff, J. Chem. Phys. 125, 194521 (2006). 68M. Khalil, N. Demirdöven, and A. Tokmakoff, J. Phys. Chem. A 107, 5258 (2003). 69T. Ishiyama and A. Morita, J. Phys. Chem. C 113, 16299 (2009). 70T. Ishiyama and A. Morita, J. Phys. Chem. C 111, 738 (2007). 45Y.

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.