Simultaneous Spectral and Temporal Analyses of Kinetic Energies in

Article pubs.acs.org/JPCA

Simultaneous Spectral and Temporal Analyses of Kinetic Energies in Nonequilibrium Systems: Theory and Application to Vibrational Relaxation of O−D Stretch Mode of HOD in Water Jonggu Jeon, Joon Hyung Lim, Seongheun Kim, Heejae Kim, and Minhaeng Cho* Department of Chemistry, Korea University, Seoul 136-701, Korea ABSTRACT: A time series of kinetic energies (KE) from classical molecular dynamics (MD) simulation contains fundamental information on system dynamics. It can also be analyzed in the frequency domain through Fourier transformation (FT) of velocity correlation functions, providing energy content of different spectral regions. By limiting the FT time span, we have previously shown that spectral resolution of KE evolution is possible in the nonequilibrium situations [Jeon and Cho, J. Chem. Phys. 2011, 135, 214504]. In this paper, we refine the method by employing the concept of instantaneous power spectra, extending it to reflect an instantaneous time-correlation of velocities with those in the future as well as with those in the past, and present a new method to obtain the instantaneous spectral density of KE (iKESD). This approach enables the simultaneous spectral and temporal resolution of KE with unlimited time precision. We discuss the formal and novel properties of the new iKESD approaches and how to optimize computational methods and determine parameters for practical applications. The method is specifically applied to the nonequilibrium MD simulation of vibrational relaxation of the OD stretch mode in a hydrated HOD molecule by employing a hybrid quantum mechanical/molecular mechanical (QM/MM) potential. We directly compare the computational results with the OD band population relaxation time profiles extracted from the IR pump−probe measurements for 5% HOD in water. The calculated iKESD yields the OD bond relaxation time scale ∼30% larger than the experimental value, and this decay is largely frequency-independent if the classical anharmonicity is accounted for. From the integrated iKESD over intra- and intermolecular bands, the major energy transfer pathways were found to involve the HOD bending mode in the subps range, then the internal modes of the solvent until 5 ps after excitation, and eventually the solvent intermolecular modes. Also, strong hydrogen-bonding of HOD is found to significantly hinder the initial intramolecular energy transfer process.

1. INTRODUCTION Kinetic energy (KE) is a fundamental dynamic observable that can be calculated at any point in time from, e.g., classical molecular dynamics (MD) simulations. While it is necessary to consider both kinetic and potential energies to describe the energetics of a given molecular system, the potential energy is often difficult to analyze because it involves many-body interactions and cannot be easily decomposed for subsets of the system. In contrast, the KE can be decomposed to atomic components and can be defined precisely for arbitrary subsets of the system. In addition, the KE is more closely associated with the dynamics of the system. For example, the classical picture of the vibrational spectroscopy is the resonance of oscillating dipoles with the electric field and the vibrational spectrum is often calculated from the linear response theory1,2 by the Fourier transformation (FT) of the time correlation function (TCF) of velocities (density of states)3,4 that can be considered as the spectral density of kinetic energy (KESD) after a mass-weighted average.5,6 In our previous paper,7 we have presented a computational method to calculate the time-dependent KESD in nonequilibrium situations. It was shown that an average KE over a given time interval is related to the FT of the velocity TCF © XXXX American Chemical Society

evaluated within the same time interval. By limiting the time span so that it does not exceed the time scale of the nonequilibrium process under study, we have obtained the nonequilibrium KESD with respect to time and frequency variables. This enabled us to calculate the time evolution of various spectral components of a molecular system and analyze the energy exchange process among those components. While this method was found to be adequate for the study of the carbonyl bond vibrational relaxation in a solvated Nmethylacetamide (NMA) molecule,7 it fundamentally relies on a short-time average and therefore its application to processes with faster time scale is expected to be problematic. In this paper, we present a new theoretical method that overcomes the limited time resolution of our previous timeaveraged KESD method. The new method employs the concept of instantaneous power spectrum introduced by Page8 and Lampard.9 Here, instead of introducing the averaging Special Issue: Jacopo Tomasi Festschrift Received: October 8, 2014 Revised: December 1, 2014

A

dx.doi.org/10.1021/jp510157y | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

time span for the TCF calculation, one of the two time points of the time correlation product of the form v(t1)·v(t2) is taken as the observation time, and the TCF is written in terms of the time variable |t1 − t2| without time average. In the original expression of Page,8 the TCF is defined such that the observation time, say t1, is larger than t2, but we also consider the opposite case where t2 > t1. These two different definitions represent the time correlation at t1 with the past and future events, respectively. In this way, an instantaneous version of the KESD can be derived, which provides the instantaneous KE at a single time t1 after integration over all frequencies. We describe formal properties of the new instantaneous KESD and optimal methods and parameters for practical applications. An independent and closely related approach has been previously reported by Yagasaki and Saito6,10 for which the present development could offer a better theoretical foundation. In contrast with the average KESD which has a limited time resolution of order 0.1 ps, the instantaneous KESD does not rely on time average and reflects a truly instantaneous system dynamics along a general nonequilibrium trajectory. Even when time average is performed on the instantaneous KESD, its spectral resolution is superior to that of the average KESD at the same level of time average. As an application of the new method, we study the vibrational relaxation of an excited OD bond of HOD in aqueous solution.11−16 In the computation, the OD bond excitation is simulated with a sudden change in atomic velocity of the oxygen and deuterium, and the velocity trajectory necessary for the KESD computation is generated by nonequilibrium MD (NEMD) simulation. We employ a hybrid quantum mechanical/molecular mechanical potential to describe the HOD and solvent H2O molecules. The QM description of the HOD molecule, implemented with the semiempirical scc-DFTB potential,17 is expected to realistically capture the vibrational anharmonicity of the HOD internal modes that is crucial for accurate simulation of the intramolecular vibrational energy redistribution process. In this study, we treat the vibrational relaxation classically, and this approach could be improved by including quantum effects18 in the future. We also perform infrared (IR) pump−probe spectroscopy experiment on the same system, i.e., 5% (v/v) HOD in water, and compare the experimental results directly with the predictions of the theory.

K av(t0 , τ ) =

=

t 0 + τ /2

∫t −τ/2

N

dt ∑ mi |u i(t ; t0 , τ )|2

0

i=1 N

∞

1 4πτ

∫−∞ dω ∑ mi |u̅ i(ω; t0 , τ)|2 i=1

(2) 1

where mi is the mass of atom i and the Parceval’s theorem is used to obtain the second equality in terms of the Fourier transform (FT) u̅i(ω;t0,τ) of atomic velocities ∞

u̅ i(ω; t0 , τ ) = =

∫−∞ dt ui(t ; t0 , τ)eiωt t 0 + τ /2

∫t −τ/2

dt vi(t )eiωt

0

(3)

We note that the summation in eq 2 may be confined to a subset of the system that we are interested in. For instance, if only the KE of solute molecules is of interest, the summation includes solute atoms. Furthermore, if necessary, we can employ a local or normal mode representation of the velocities instead of the atomic ones with corresponding changes in the mass mi. From eq 2, we can identify the spectral density Φav(ω;t0,τ) of KE (KESD) Φav (ω; t0 , τ ) =

1 4πτ

N

∑ mi |u̅ i(ω; t0 , τ)|2 i=1

(4)

which provides the average KE Kav(t0, τ) when integrated over the entire frequencies. This equation guarantees Φav(ω;t0,τ) to be non-negative. To proceed further, we introduce the TCF Ci(t;t0,τ) of velocities Ci(t ; t0 , τ ) = ⟨u i(T + t ; t0 , τ ) ·u i(T ; t0 , τ )⟩ (t0 − τ /2 ≤ T ≤ t0 + τ /2 − t )

(5)

and assume that its T-dependence is negligibly small. This is the key approximation in our previous theoretical development of nonequilibrium KESD.7 It requires that the time span τ be shorter than the time scale of the nonequilibrium process under study. Then, the KESD in eq 4 can be expressed in terms of Ci(t;t0,τ) after ensemble average Φav (ω; t0 , τ ) =

1 4π eiωt

2. THEORY AND COMPUTATION

N

∑ mi ∫ i=1

∞

−∞

dt(1 − |t | /τ )Ci(t ; t0 , τ ) (6)

In our previous paper, we ignored the term |t|/τ in the integrand, which is justified when ∑Ni=1miCi(t;t0,τ) decays fast compared to τ.7 Equation 6 shows that the KESD averaged over the time span τ can be calculated from the FT of velocity TCF after applying a triangular window function signified by the factor (1 − |t|/τ), which is a direct consequence of the limited time range defined by eq 2. We note that the presence of this factor does not affect the integrated KESD over all frequencies since the integral depends only on the value of the integrand at time zero. It also means that the total average KE is not affected by the choice of the window function often used in the FT process19 as long as it has a value of one at time zero. Hereafter, we will refer to Φav(ω;t0,τ) as “aKESD”. We note that shorter time span τ is desirable for better time resolution of aKESD and also for the validity of the approximation made in eq 5.

2.1. Average Spectral Density of Kinetic Energies (aKESD). We first present a brief description of the theory of average KESD in nonequilibrium situations.7 Let us consider a time series of atomic velocities ui(t;t0,τ) defined within a time span τ centered at t0 ⎧ v (t ) t0 − τ /2 ≤ t ≤ t0 + τ /2 u i (t ; t 0 , τ ) = ⎨ i ⎩ 0 otherwise

1 2τ

(1)

where vi(t) is the atomic velocities, e.g., from a general nonequilibrium MD (NEMD) trajectory. The KE of the system, centered at time t0 and averaged over the time span τ, is then given as B



Article

over a finite time window and the result is defined at a single time point of T. Although the preceding derivation and the original development of Page8 consider only the correlation of the present-time variable with those in the past, it is also possible to derive an analogous expression in terms of the correlation with the future. In this case, we begin with the time span T ≤ t ≤ t0 and corresponding velocities,

However, it cannot be made too small for a reliable FT in eq 6 and it is usually necessary to find an optimal value of τ on a system-by-system basis. In equilibrium situations, we can take the limit τ → ∞ in eq 1 and Ci(t;t0,τ) becomes independent of t0, τ, and T in eq 5. Then, the KESD in eq 6 becomes Φ′eq (ω) =

N

1 4π

∑ mi ∫

∞

−∞

i=1

dtCi(t )eiωt

⎧ v (t ) T ≤ t ≤ t 0 u i (t ; t 0 , T ) = ⎨ i ⎩ 0 otherwise

(7)

7

which can also be derived directly using the Wiener−Khinchin theorem.1,2 2.2. Instantaneous KESD: A New Approach. To overcome the limitations of the aKESD approach, we have developed a new approach to a truly instantaneous KESD with unlimited time resolution. This is based on the notion of instantaneous power spectrum introduced by Page8 and Lampard.9 We first define a time span beginning at an arbitrary time t0 in the past and ending at time T at which we would like to find the KESD. The atomic velocities in this time span are expressed as ⎧ v (t ) t 0 ≤ t ≤ T u i (t ; t 0 , T ) = ⎨ i ⎩ 0 otherwise

Then, the instantaneous KESD at time T can be obtained as t0− T mi dτ vi(T ) ·vi(T + τ ) cos ωτ 2π 0 ∞ m d τ u i (T ; t 0 , T ) · = i 2π −∞ u i(T + τ ; t0 , T ) cos ωτ

ρi F (ω; T , t0) =

(14)

When integrated over all frequencies, ρFi (ω;T,t0) gives the same instantaneous KE at time T shown in eq 12. The superscript ‘F’ in ρFi (ω;T,t0) denotes that the velocity correlation at time T is calculated using the f uture values at T + τ. Hereafter, ρFi (ω;T,t0) will be termed “fKESD”, and the name “iKESD” will be used to collectively refer to the two instantaneous KESDs, pKESD and fKESD. In terms of the non-negative wavenumber ν̅ instead of the positive/negative angular frequency ω, the two iKESDs can be written as

(8)

∞

=

∫−∞ dt ui(t ; t0 , T )eiωt ∫t

T

dt vi(t )eiωt

(9)

0

ρi F/P (ν ̅ ; T , t0)

We also introduce the power spectrum of velocities Pi(ω;t0,T) defined as Pi(ω; t0 , T ) =

∫0

= 2cmi Re[

1 |u i(ω; t0 , T )|2 2π ̅

1 = 2π

∫t

T

dt1 vi(t1) · 0

t1− t 0

∫t −T

dτ vi(t1 − τ )e

∫0

iωτ

mi ∂ Pi(ω; t0 , T ) 2 ∂T T − t0 m dτ vi(T ) ·vi(T − τ ) cos ωτ = i 2π 0 ∞ m d τ u i (T ; t 0 , T ) · = i 2π −∞ u i(T − τ ; t0 , T ) cos ωτ

ρiP (ω; T , t0) =

∫ ∫

(11)

ρPi (ω;T,t0)

since integrating over all frequencies yields the instantaneous KE of atom i at time T mi |vi(T )|2 2

d τ u i (T ; t 0 , T ) ·

∞

F/P dνρ ̅ i (ν ̅ ; T , t0) =

mi |vi(T )|2 2

(15)

where c is the speed of light and the superscripts “F/P” and the positive/negative sign in front of τ corresponds to fKESD/ pKESD, respectively. Similarly, the aKESD in eq 6 can be redefined in terms of non-negative ν̅.7 In equilibrium, we can freely choose the time span for the KESD evaluation. Thus, we carry out the time average by integrating ρF/P i (ν;T,t ̅ 0) over T between 0 and T′, apply the ensemble average to the time correlation product ui(T;t0,T)·ui(T ± τ;t0,T), and take the limit t0 → ± ∞ and T′ → ∓∞ for the fKESD and pKESD, respectively. Then, it can be shown that the time and ensemble averaged fKESD and pKESD are reduced to the equilibrium KESD in eq 7. We note that the present derivation of instantaneous KESDs does not involve any ensemble average unlike the aKESD case. However, the ensemble average can greatly improve computational reliability since the velocity correlation product of the form ui(T;t0,T)·ui(T ± τ;t0,T) converges to zero much faster with it. In this case, the calculated KESDs become nearly independent of the maximum correlation time |T − t0| as long as it is larger than the decay time scale or correlation time of the velocity TCF. We also note that these instantaneous KESDs may take negative values in some frequency regions, due to the lack of the time averaging in the underlying KE shown in eq 15. This is to be contrasted with the aKESD, whose non-negative

1

The instantaneous KESD at time T is given by the time derivative of Pi(ω;t0,T)

∞

∞

u i(T ± τ ; t0 , T ) cos ωτ ]

(10)

∫−∞ dωρiP (ω; T , t0) =

∫ ∫

and their FT, ui̅ (ω;t0,T), are given as u̅ i(ω; t0 , T ) =

(13)

(12)

ρPi (ω;T,t0)

The superscript P in denotes that the velocity correlation at time T is calculated using the velocities in the past, T − τ, in eq 11. Hereafter, this KESD will be referred to as “pKESD”. Note that we completely avoid time-averaging of KE C



Article

value better represents the excitation energy absorbed by an infrared photon in actual experiments. However, the resulting KESD with the latter showed significant red shift of the OD band due to classical anharmonicity effect. Such undesired behavior can be effectively reduced by using smaller excitation energy. Other than this effect related to the classical anharmonicity, we have not observed any significant qualitative difference in the relaxation behavior from the two different perturbations as will be shown in the Results section. This indicates that the present vibrational relaxation process conforms to linear response theory1,2 over a wide range of perturbation magnitude. The simulation system contains one HOD molecule described by the scc-DFTB QM potential17 and 277 H2O molecules described by the flexible SPC/Fw MM potential.21 The system is under periodic boundary conditions. To assess the reliability of the employed potentials, peak positions of the vibrational spectrum of the hydrated HOD molecule were first determined by the FT of the HOD dipole TCF (Table 1).1 It

character stems from the time averaging in the KE indicated in eq 2. While the piecewise negative spectral density is mathematically correct, it causes some practical difficulties in its interpretation. In test calculations, we observe negative regions in iKESD routinely from a single trajectory, even in equilibrium. If we ensemble/time average the iKESD, however, the negative regions tend to disappear in most cases (see below). In addition, the fKESD and pKESD result in the same instantaneous KE when integrated over all frequencies, but they are generally not identical by themselves. We have found from the application to hydrated HOD system described below that their difference is largest right after an abrupt perturbation to the OD bond velocity, and it gradually decreases as the system approaches equilibrium. This indicates that the oppositely directed time correlations in fKESD and pKESD may be related to the extent of nonequilibrium nature of the system. This nonequilibrium dynamic aspect of the time correlations and instantaneous spectral densities merits further formal investigation in the future. 2.3. Computational Details. We calculated the KESDs described above for the vibrational energy relaxation process of the excited OD bond of HOD in water. The process was simulated with the NEMD method in which the OD bond is excited at time zero by sudden changes in the atomic velocities. Specifically, we assume that (i) only the velocities of the oxygen and deuterium are changed initially due to bond excitation such that the bond KE increases by ℏωvib, (ii) the perturbation is along the original OD bond vector direction, and (iii) the center of mass of the bond and its velocity are unchanged by the perturbation. Under these conditions, the perturbed velocities v are given in terms of the unperturbed velocities v0 by vD = v 0D + εDr0̂ ,

Table 1. Peak Frequencies (in cm−1) of the Intramolecular Modes in the Vibrational Spectra of the HOD System Hydrated by H2O Moleculesa system QM HOD/MM H2O MM HOD/MM H2O NMM HOD/ NMM H2O QM HOD/ NMM H2O exp. ([D2O] = [H2O])b

0 vO = v O + εOr0̂

T (K)

HOD bend

OD stretch

OH stretch

298

1368 (+41)

2649 (−151)

3649 (−229)

297

1295 (+90)

2648 (−8)

3634 (−16)

297

1368 (+85)

2644 (−13)

3634 (−16)

299

1374 (+47)

2658 (−142)

3645 (−233)

1450 (+47)

2495 (−229)

3415 (−293)

a

(16)

Results for a few combinations of the QM and MM potentials are shown, and the peak shifts from the corresponding isolated HOD molecule are indicated in parentheses. The QM HOD/NMM H2O combination was used in the production simulation. In all cases, the reduced LJ radius for the QM oxygen was used as described in text. (QM: scc-DFTB potential, MM: SPC/Fw potential, NMM: SPC/Fw with increased bending force constant). bReference 29.

where r̂0 is the unit vector parallel with rD − rO and μ is the reduced mass of the OD bond. This method of perturbation is different from our previous approach7,20 to carbonyl bond excitation in that only the bond velocity and not the atomic positions are changed. This difference is not expected to significantly affect the computational results because, in the present classical approach, the initially deposited vibrational kinetic energy will be quickly repartitioned between the kinetic and potential energies of the local oscillator within a few cycles of vibration (tens of fs) and this process is effectively decoupled from the major relaxation mechanisms which have 100 fs or longer time scales (see below). In addition, the present velocity perturbation scheme has the advantages that the excitation energy uniquely determines the changes in velocities without the use of random numbers and the (harmonic) functional form of the bond stretch force field does not need to be assumed. We have performed two sets of NEMD simulations with OD bond excitation energies ℏωvib of 0.592 and 7.660 kcal/mol, corresponding to frequencies of 207 and 2679 cm−1. The former is the thermal energy at 298 K and the latter corresponds to the peak frequency of the OD stretch band determined from the FT of the HOD dipole TCF. The latter

showed reasonable agreement with experimental peak positions, but the scc-DFTB potential somewhat underestimates the HOD bending mode frequency and overestimates the OD and OH stretch frequencies. Furthermore, another system with an HOD described by the MM potential instead of the QM one showed that the OD and OH stretch frequencies of the two models coincidently agree within 15 cm−1 but their HOD bending frequencies differ by ∼70 cm−1. This finding indicates that the internal consistency of the present QM/MM description can be greatly improved by a simple modification to the bending potential of the MM model. Therefore, we have increased the bending force constant of the SPC/Fw model by 13% to 86.00 kcal/(mol·rad)2 so that two HOD molecules described respectively by the QM and the modified MM potentials (NMM in Table 1) exhibit nearly the same frequencies for all the internal modes. Note that the modified MM potential is applied to solvent H2O molecules, and no change is made to the scc-DFTB potential used for the QM HOD molecule. In addition, following our previous study,22 we have used a new Lennard-Jones (LJ) radius for the QM oxygen, which is 6% smaller than the standard value of the SPC/Fw

εD = (p − p0 )/mD , p=

p0 |p0 |

εO = −(p − p0 )/mO

(p02 + 2μℏωvib)1/2 ,

0 ) ·r0̂ p0 = μ(v 0D − v O

D



Article

bond can reveal spectral compositions of the mode which could otherwise be obtained only through normal-mode analysis. In principle, the iKESD can be calculated at any T after the OD bond excitation by choosing the other end point t0 of the time span in eqs 8 or 13 such that the velocity TCF converges well to zero within |T − t0|, that is, |T − t0| should not be too small compared to the self-correlation time of velocities. Since we choose to take time averages of the iKESD over 100 fs of T as mentioned above, we can obtain the iKESD beginning at ∼50 fs after the bond excitation is applied. However, we have found that the pKESD exhibits unphysical features if the times prior to the bond excitation are included in the TCF evaluation, due to the discontinuity in the trajectory introduced by the perturbation. Thus, we have decided to carry out the analyses only with the fKESD that depends only on the postexcitation trajectories and does not suffer from this artifact. At times when both pKESD and fKESD can be evaluated reliably, their difference is not significant enough to change the interpretation of the result, justifying the use of fKESD at all times. The optimal value of |T − t0| for all values of T is ∼1 ps, but it was found that |T − t0| of 400 fs can also be used without introducing too much change in the band envelops if the Kaiser window function19 with θ = π is applied to the TCF prior to the FT. We used the latter protocol in the analyses presented below. The FT was performed with the maximum entropy method24 since it provides sharper signals with smaller side lobes than with the conventional discrete FT.

model. This strengthens the QM−MM interaction to the level of the MM−MM one and improves the consistency in intermolecular interactions. Vibrational peak positions of the various combinations of the potential models are summarized in Table 1. The system was first equilibrated under constant NPT and constant NVT conditions for 9 ns in total, using the time step size of 0.5 fs, the particle mesh Ewald method for long-range electrostatic interactions, and an LJ cutoff distance of 9 Å. The density of the system was found to be 0.995 g/cm3 at the end of this equilibration. Then, the canonical ensemble of the system was prepared from a constant NVT simulation with Langevin thermostat at 298 K and 500 or 1000 samples from this run were used as initial configurations for the NEMD simulation. In the NEMD stage, each of the 500 or 1000 samples was velocityperturbed according to eq 16 and then propagated for 20 ps under the constant NVE condition. The nonequilibrium trajectories were saved in 4 fs interval for analysis. The Amber 11 package was used in all MD simulations.23 The average and instantaneous KESDs were calculated from the NEMD velocity trajectories. Test calculations have confirmed that all integrated KESDs perfectly reproduce the instantaneous or short-time averaged KE as required by eqs 2 and 12. It was also found that the iKESDs evaluated in 4 fs intervals of T from a single trajectory can reliably capture individual KE variations of the OD bond within a single cycle of vibration of ∼13 fs. This is in stark contrast with the aKESD for which the averaging time span τ of at least 100 fs is necessary to obtain a reliable signal. While this fine time resolution of the iKESD has a great potential as an analytical method, it provides too much detail for the present purposes. Thus, in the following application, we have taken a time average of iKESD with respect to T over a time span of 100 fs. This choice provides a sufficient average of the OD stretch mode, which is the main focus of this study. In addition, ensemble averages were taken of the underlying velocity TCF over the 500 or 1000 samples. We also note that the iKESD in general produces much sharper signal than the aKESD at the same level of time average. We will therefore use the iKESD exclusively in the following analyses. Since the bond excitation scheme employed here affects only the OD bond velocity of the HOD molecule, KESDs for the local OD stretch mode provide the most direct measure of the vibrational energy relaxation. From the OD bond distance rb = | rD − rO|, the bond velocity vb is given by v b = rḃ = (vD − vO) ·r0̂

3. EXPERIMENT The sample of 5% (v/v) HOD solution was prepared by mixing 2.5% D2O molecules purchased from Sigma-Aldrich (99.9 atom %) with pure water obtained from the Milli-Q water purification system. The infrared (IR) pump−probe (PP) measurements were performed using mid-IR pulses.25 Although essentially the same experimental studies were performed and reported before by other groups,11−13 we have carried out the present IR PP measurements for direct and quantitative comparisons with the present nonequilibrium MD simulation results on the vibrational energy relaxation of the IR-pumped OD stretch mode of HOD in water at room temperature. Here, the 800 nm beam was generated by Ti:sapphire oscillator/ regenerative amplifier (Tsunami, Spitfire; Spectra Physics) and the beam pumped optical parametric amplifier to produce two near-IR pulses centered at ∼7500 and ∼5000 cm−1. The resulting two pulses were difference-frequency-mixed in AgGaS2 so that mid-IR pulses with center frequency of 2530 cm−1 were produced, where pulse duration is ∼70 fs and pulse energy is ∼1 μJ. In the isotropic signal from the experiment, a persistent spectrum is observed at long delay time. The origin of this nondecaying component is identified as the temperature rise in the sample caused by the absorption and subsequent thermalization of the energy of the pump pulse.11,12 This heating effect results in the thermalization component that contributes to the observed transient spectra of HOD together with the pure decay component. The two components were determined from the total signal by fitting to the three-state kinetic model of Bakker and co-workers,12,13 which presupposes that the transition rate constants between the states are frequency-independent

(17)

The OD bond KESDs were calculated using only vb in eqs 5, 6, and 15 with necessary changes such as changing vector inner product to scalar product and using the reduced mass μ instead of the atomic mass. We have also calculated KESDs for the HOD molecule using the three atomic velocities, and the total KESDs including all atomic velocities of the system. They are respectively used to investigate the intramolecular vibrational energy redistribution and the intermolecular vibrational relaxation pathways. We note that the total KESD contains information on the entire system dynamics including those conveyed by OD bond/HOD KESDs. However, just as in experimental spectroscopy, the total KESD may suffer from signal overlaps and congestions, and it is sometimes advantageous to separately analyze the KESDs for subsets of the system, i.e., OD group, HOD molecule, and solvent water molecules. In addition, the local mode KESD such as the OD E



Article

⎡ ⎤ τ τ1 * e−t / τ* − Δα(ω , t ) = Δα(ω , ∞)⎢1 + e−t / τ1⎥ τ1 − τ τ1 − τ ⎣ ⎦ * * + N1(0)[σ12(ω) − 2σ01(ω , 0)]e−t / τ1

attributed to the new equilibrium state generated by the dissipation of the coherent vibrational energy absorbed from the pump pulse. In fact, Fayer and co-workers11 have shown that the long time spectrum is very similar to the temperature difference spectrum corresponding to 2 K increase in the system temperature, supporting its thermal origin. Following the works of Fayer and Bakker, we have attempted a separation of the thermalization component by employing the three-state kinetic model of the Bakker group introduced in the previous section. The isotropic IR PP signal was globally fit to the kinetic model expression in eq 18 using the entire experimental data both in the time and frequency domains. From this, we can obtain the thermalization and pure decay components, corresponding to the first and second terms of eq 18, respectively. These two components are displayed in Figure 1b,c. At small T, the pure decay component in Figure 1b closely resembles the total isotropic signal in Figure 1a. However, the decay pattern of the former is uniform over the probe frequencies due to the assumptions of the model. The thermalization component in Figure 1c grows from zero and approaches the long time isotropic signal at large T. From the fitting procedure, the decay time τ1 of the pure decay component is found to be 1.794 ps, and the decay time τ* of the thermalization component, 0.825 ps. These values agree well with those reported by Rezus and Bakker.12 In Figure 2a, the time profile of the isotropic IR PP signals at several different frequencies are plotted together with the fitting curves from the kinetic model. For better comparison, the data were normalized to the values at the earliest time point of 0.2 ps. The excellent

(18)

where Δα(ω,t) is the isotropic PP signal at time t, the subscripts 1 and * denote the pure decay and thermalization components respectively, N1(0) is the initial population of the singly excited state, σ01(ω,t) and σ12(ω) are the cross sections of the 1 ← 0 and 2 ← 1 transitions, respectively, with the latter assumed to be time-independent. The first term in (18) is the thermalization component that gives the persistent PP signal in the long time limit and the second term is the pure decay component. Note that the thermalization component has two time constants with opposite signs whereas the pure decay term is described by a single exponential decay. We emphasize that the time-dependence of each component is identical regardless of the signal frequency.

4. RESULTS 4.1. Experimental Results. The isotropic signal from the polarization-controlled IR PP experiment is presented in Figure 1a. The signal is mainly composed of the positive 1 ← 0 band

Figure 1. (a) Isotropic signal from the vibrational pump−probe experiment on the hydrated HOD system. (b) Pure decay component of the signal obtained from the three-state kinetic model. (c) Thermalization component of the signal from the same model.

Figure 2. (a) Time profile of the isotropic pump−probe signal at several frequency points in the 1 ← 0 transition region (empty symbols) of the OD stretch band. The signals are normalized to the values at 0.2 ps. Solid lines are the result of fitting to the kinetic model. (b) Normalized time dependence of the pure decay and thermalization components determined from the fitting. The time dependence of the experimental isotropic PP signal is well described by these two kinetic components at all frequencies with different linear combinations. The pure decay component has a single exponential decay time of 1.794 ps. The thermalization component contains a rise term with the same 1.794 ps time constant and a decay term with a 0.825 ps time constant. See text for details.

of the OD stretch mode above 2430 cm−1, with a smaller negative band corresponding to the 2 ← 1 transition around 2400 cm−1. The main signal decays with T reflecting the population loss of the singly excited state generated by the pump pulse. The decay behavior is not uniform over the frequency range shown. In addition, the spectral shape of the IR PP signal after ∼8 ps appears to be stationary without further decay. This persistent spectrum at large T was previously F



Article

agreement of the two indicates the validity of the kinetic model. The normalized time profiles of the two kinetic components are shown in Figure 2b. With a suitable frequency-dependent linear combination according to eq 18, these two components can completely describe the apparent kinetics of the isotropic IR PP signal. We have also tested a different approach to analyze the time dependence of the PP signal by attempting a multiexponential fitting of the signal at a few individual frequencies. This fitting resulted in a common major decay component of ∼1.8 ps time scale at many of the frequencies but smaller components were not consistent over the frequencies. This finding casts some doubt on the robustness of the preceding result. We also considered the total signal by first integrating the isotropic PP spectrum over the 1 ← 0 band between 2450 and 2650 cm−1 and then fitting the integrated data to a threeterm exponential function. This procedure resulted in two decay components with time constants of 0.862 ps (14.3%) and 1.822 ps (85.7%), in good agreement with the decay times of thermalization and pure decay components obtained above. Therefore, as expected from the assumptions of the model, the kinetic model analysis presented above reflects only the average kinetic behavior of the spectral band. These results suggest that, even though the global fitting of the present type agrees very well with the time- and frequency-resolved PP signal, the underlying kinetic model needs to be interpreted more carefully. We note that the analysis presented above due to Rezus and Bakker12 is similar in spirit to the approach by Steinel, Fayer, and co-workers11 in that both of them assume f requencyindependent kinetics of the underlying components and introduce two types of kinetic constants. However, the model of Bakker et al. assumes the presence of the intermediate state, and the decay kinetics of the pure decay component also contributes to the rise of the thermalization component in an important way. On the other hand, the approach by the Fayer group is based on the singular value decomposition analysis, and the decay and rise kinetics are apparently not related to each other. 4.2. Computational Results. In the following, we present the computational analysis of the NEMD simulations on the OD stretch relaxation process. We first study the overall features of the calculated KESD in equilibrium. Figure 3a displays the equilibrium iKESD for the OD bond, HOD molecule, and the entire system calculated from the equilibrium trajectories after time and ensemble averages. The OD bond iKESD shows a single dominant peak at ∼2700 cm−1 arising from the OD stretch degree of freedom. We can also identify much smaller peaks at 3765 and 1370 cm−1, corresponding to OH stretch and HOD bending motions, which arise from small discrepancy between the OD local mode and the three normal modes observed. We have obtained the KE of each band by integrating the iKESD over the frequency ranges specified in Table 2. From the entries for OD bond in Table 2, the contribution of the 2700 cm−1 band to the total is 96.6% (0.283/0.293). We emphasize that this kind of analysis can be performed without the explicit use of the potential energy information unlike the normal-mode analysis. The baseline of the OD bond iKESD in Figure 3a is ∼1 × 10−6 kcal·cm/mol, which is the limit imposed by the “noise” parameter used in the maximum entropy method. The iKESD for the HOD molecule in the same figure shows three major peaks already observed in the OD bond iKESD and also a broad band below 1200 cm−1 arising from the hindered translation and rotation of the HOD

Figure 3. (a) Time-averaged instantaneous KESD of the hydrated HOD system in equilibrium. KESDs for three subsets of the system are shown. (b) Nonequilibrium fKESD of the HOD molecule as a function of the elapsed time T since the OD bond excitation by 7.660 kcal/mol. fKESD at each T is the result of a time average over a 0.1 ps segment centered at that T.

Table 2. Band Kinetic Energies from Integration of the iKESD in Equilibrium kinetic energy (kcal/mol) −1

band (cm ) 0−360 360−1200 0−1200 1200−1475 1200−1750 2370−2880 3400−4000 0−4170

(HOD bend) (H2O/HOD bend) (OD stretch) (OH stretch)

OD bond

HOD

system

0.283 0.293

1.798 0.253 0.279 0.282 2.665

268.370 226.642 495.012 79.578 0.934 158.756 744.112

molecule as a whole. We note that the OD bond iKESD is slightly larger than the HOD iKESD near the blue side of the OD stretch band in the 2700−2750 cm−1 region. This anomaly is also reflected in the integrated iKESD in Table 2, which shows the OD stretch band KE of 0.283 kcal/mol for the OD bond iKESD and 0.279 kcal/mol for the HOD iKESD. We ascribe this small discrepancy to computational error, which does not significantly affect the interpretation of the results. Finally, the iKESD for the entire system in Figure 3a shows a strong HOH bending band at ∼1550 cm−1 and an OH stretch band at ∼3630 cm−1 that reflect the internal motion of the MM H2O molecules. From the band KE in Table 2, we note that the energy content of the low frequency band below 1200 cm−1 is about 2/3 of the total, both in HOD and the entire system, reflecting the fraction (6/9) of the translational and rotational degrees of freedom of a single water molecule. An example of the time-dependent iKESD in nonequilibrium situations is presented in Figure 3b, which shows the fKESD of the HOD molecule calculated from the NEMD trajectories for the large perturbation case after a time average over a 0.1 ps range and an ensemble average over 500 independent NEMD trajectories. The figure captures the gradual decrease of the OD G



Article

for the solvent KE shows that the energy transfer to solvent is rather slow until ∼1.5 ps but it increases rapidly in the 2−3 ps range before converging to a new equilibrium about 2.5 kcal/ mol higher in energy than the pre-excitation state. An analogous analysis for the small perturbation case results in two decay components of the OD bond KE with time constants of 0.08 ps (13.6%) and 2.27 ps (86.4%) that are very similar to the values found above for the large perturbation case. This indicates that the decay kinetics of the OD bond KE is not significantly affected by different perturbation magnitudes. In the small perturbation case, however, the time profiles of the KEs for the rest of the HOD molecule and the solvent are severely masked by thermal fluctuations, and their time constants are difficult to extract reliably. 4.3. Simultaneous Spectral and Temporal Analysis. We now turn to a more detailed analysis of the iKESD in both the time and frequency domains. We first focus on the relaxation of the OD bond KE using the OD bond fKESD. In Figure 5, the OD bond fKESDs after a large (7.660 kcal/mol)

stretch band after its excitation at time zero and also the intramolecular vibrational energy redistribution (IVR) processes of smaller magnitude taking place mainly in the bands of the HOD bending and OH stretch motions (see below). The nonequilibrium fKESD after 10 ps closely matches the equilibrium iKESD, indicating that the vibrational energy transfer from the HOD molecule to solvent is almost complete within 10 ps. The time-dependent changes of the KE of the OD bond, the HOD molecule, and the solvent H2O molecules are displayed in Figure 4 for the large perturbation case. They were calculated

Figure 4. Changes in kinetic energies of the OD bond (OD-s), the rest of the HOD molecule (HOD excl. OD-s), and 277 solvent H2O molecules (Solvent) after the OD bond excitation by 7.660 kcal/mol at time zero. The solid lines were obtained directly from the atomic kinetic energies with 0.1 ps moving averages and the data points marked with empty symbols represent the integrated fKESD with time average over 0.1 ps segment. The filled symbols at time zero represent the initial values in the absence of the perturbation. In all cases, the average equilibrium values in Table 2 were subtracted from the data.

from the atomic velocities, and the curve labeled ‘HOD excl. OD-s’ was obtained by subtracting the OD bond KE from the HOD one at each time (solid lines). In the same plot, we have also included the KEs obtained by integrating the nonequilibrium fKESDs over the entire frequency range from 0 to 4170 cm−1 (empty symbols). The KEs from these two different methods perfectly agree with each other. The average equilibrium values were subtracted from each data, and the curves starts at 48 fs due to the lead time necessary for a 0.1 ps time average performed. The starting value of the OD stretch (OD-s) KE is 3.654 kcal/mol, about one-half of the initially deposited energy. The other half is accounted for by the increase in the bond potential energy during the 0.1 ps averaging time span that we do not monitor directly. The decay times of the OD-s KE were determined as 0.12 ps (9.6%) and 2.42 ps (89.9%) from a three-term exponential fit. The KE for the rest of the HOD molecule (HOD excl. OD-s) directly reflects the IVR process, and it shows the typical rise-decay pattern. A three-term exponential fit of this curve resulted in two rise components with time constants of 0.14 ps (amplitude −0.871 kcal/mol) and 0.65 ps (−0.378 kcal/mol) and a single decay component of 2.61 ps (amplitude 1.32 kcal/mol). The first rise time scale of 0.14 ps matches the fast component (0.12 ps) of the OD bond KE decay. The two rise components indicate that the IVR process takes place well within subpicosecond time scale. The intermolecular vibrational energy relaxation (VER) of the HOD molecule, reflected in the 2.61 ps decay component, takes only slightly longer than the relaxation of the OD bond KE (2.42 ps). Finally, the curve

Figure 5. Relaxation of the OD bond fKESD after bond excitation by (a) 7.660 kcal/mol and (b) 0.592 kcal/mol. In both cases, the equilibrium KESD was subtracted from each fKESD.

and a small (0.592 kcal.mol) perturbation are presented at multiple time points. For clarity, we have subtracted the equilibrium iKESD from each curve. The resulting ΔKESD shows gradual decrease with time, but its frequency dependence is rather nonuniform. In particular, in Figure 5a, we note the decay up to 0.5 ps proceeds slightly faster on the blue side of the band, as can be seen from the asymmetric line shapes at 0.5 and 1.0 ps. Between 1.0 and 1.5 ps, however, a kind of transition in mechanism takes place, and the amplitude in the red side begins to decrease faster than in the blue side, resulting in the crossing of the curves at 1.0 and 1.5 ps. This kind of uneven decay continues until the fKESD reaches equilibrium at ∼10 ps. Comparing the initial fKESD at 48 fs and the later ones between 1.5 and 5 ps, we can clearly identify a blue shift of the band. In the small perturbation case shown in Figure 5b, the early decay pattern is analogous to the large perturbation case, and we can find the same doublet feature at 0.7 ps as in Figure 5a at 1.0 ps. However, the later part of the decay does not show pronounced red shift as was found in the large perturbation H



Article

the experimental finding that the vibrational relaxation rate does not strongly depend on probe frequency. Nonetheless, the analysis presented above also demonstrates that the apparent frequency-dependent decay of the OD bond KE is rather complicated. At short times, the decay is faster on the blue side of the band, which is associated with fewer hydrogen bonds (H-bonds) around the chromophore.26,27 This indicates that the stronger H-bonding of the low-frequency ensemble interferes with efficient IVR process. On the other hand, the long time decay beyond ∼1.5 ps is found to be faster in the red side of the band. From the comparison of the large and small perturbation cases, we can unambiguously attribute this to the classical anharmonicity manifested by the excited OD bonds, which is not expected to be observable in a realistic quantum mechanical system. If this effect is minimized by reducing the excitation energy (small perturbation case), we obtain more or less frequency-independent decay of the signal, fully consistent with existing experimental studies. We note that the frequency-dependence of the type discussed here is determined by the relative rates of the vibrational relaxation and the spectral diffusion/H-bonding dynamics,11 and it is manifested differently, e.g., in the OH bond vibrational relaxation in H2O or D2O solvents.14,28 To study the effect of the initial H-bonding environment on the vibrational relaxation of the OD bond, we have classified the entire nonequilibrium ensemble according to the number of H-bonds formed by the HOD molecule at time zero. We determined the H-bond number using the criteria that the O− O distance is smaller than 3.4 Å and the O−H(D)-O angle is larger than 135 degree. In Figure 7, the OD bond KEs of the

case. We can attribute this red shift of decay to the anharmonicity of the OD stretch mode described by the sccDFTB potential. As a result of the initial excitation, the system is placed in high energy states, which are much more anharmonic than near the potential minimum, leading to the red shift of the initial fKESD. As the system relaxes to lower energy states that are more harmonic, the frequency of the band would shift to the blue side. This also explains the smaller degree of blue shift observed in the small excitation case. We emphasize that this kind of anharmonicity effect is purely classical in nature. In quantum systems, the energy states are discrete, and the effect of anharmonicity is hardly manifested unless the 2 ← 1 transition is monitored together with the 1 ← 0 transition. In this regard, the small perturbation case in Figure 5b is believed to be more realistic, even though its excitation energy is much smaller than in experiments. To further investigate the apparently frequency-dependent decay pattern observed in the fKESD, we have plotted in Figure 6 the time profiles of the fKESD along several frequencies. In

Figure 6. Time profile of the difference KESD shown in Figure 3 at a few frequencies. The results for bond excitation energy of (a) 7.660 kcal/mol and (b) 0.592 kcal/mol are shown. The data are normalized to the values at T = 48 fs.

Figure 7. Decay of the OD bond kinetic energy of subensembles classified by the initial number of H-bonds (NHB) formed by the HOD molecule. The energy of bond excitation at time zero was 7.660 kcal/ mol. Each curve is the result of 0.1 ps moving average.

the large perturbation case of Figure 6a, the initial decay up to 2 ps is indeed faster in the low frequency region (disregarding the 2680 cm−1 curve, which is affected more from thermal fluctuations due to its smaller unnormalized amplitude), but the later decay is faster in the high frequency region with the crossover between 1.5 and 2.0 ps, as was found in Figure 5a. In the smaller perturbation case shown in Figure 6b, the decay curves appear to be identical if thermal fluctuation is accounted for. We have attempted a three-term exponential fitting of the decay curves in Figure 6a. Due to large fluctuations observed at short times, the short time components were not very reliable, but the long time component showed nearly monotonic increase from ∼1.2 ps at 2480 cm−1 to ∼3.0 ps at 2640 cm−1, indicating more than 2-fold increase in the relaxation time of the blue side signal. On the other hand, the curves for the small perturbation case in Figure 6b yielded much smaller variation of decay times in the 2.0−2.7 ps range, which is consistent with

subensembles with different H-bond numbers are displayed as a function of time. The figure clearly shows that stronger Hbonding leads to slower decay of the OD bond KE. This is qualitatively consistent with the short-time behavior of the fKESDs observed in Figure 6a, where the red sideband of the spectrum is found to decay more slowly than the blue side. However, we emphasize that Figure 7 displays the bond kinetic energy rather than the spectral density at a given frequency as in Figure 6. Because of this difference, the different initial decay behavior of the subensembles persists until the equilibrium is reached without the crossover observed in Figure 6a. 4.4. Energy Transfer Pathways. Additional important information conveyed by the iKESD is the energy contents of different spectral regions and their time evolution. Using the three types of the nonequilibrium fKESDs for the OD bond, I



Article

(68% of the initial OD bond KE). We can see from this figure that more energy is transferred to the low frequency band below 360 cm−1 than to the libration band in the 360−1200 cm−1 region. Collectively, the energy transfer pathways can be deduced from the band KE variations in time, and the result is shown in Figure 9. It shows that (i) ∼20% of the excess KE of

the HOD molecule, and the entire system, we have calculated the band KEs over the frequency ranges specified in Table 2. By subtracting the corresponding equilibrium band KEs also shown in Table 2, and taking the differences among them, we have calculated the time evolution of the KEs of various subsets of the system. The result of this analysis for the large perturbation case is displayed in Figure 8 for bands with

Figure 9. Vibrational energy transfer pathways of the HOD in water deduced from the band KE evolution. On the left side, the band frequency range is indicated in cm−1 units. The values in the boxes are the percentage of the transferred kinetic energy relative to the excess KE of the OD bond at 48 fs, 3.654 kcal/mol. Numbers in the shaded boxes represent the energy loss of the band, and those in empty boxes, the energy transferred to the band. The solid arrows represent the energy flow from the OD stretch band, and the dashed arrows, the secondary flow originating from the rest of the system. For clarity, only pathways involving more than 3% of the total excess energy are shown.

the OD bond is initially transferred to the HOD bending mode within 1 ps, (ii) the intramolecular modes of the solvent H2O (OH stertch and HOH bending modes) transiently take up more than one-quarter of the energy either directly from the OD bond or from the HOD bending mode over the 0−5 ps range, (iii) eventually more than 75% of the excess KE is stored in the low frequency solvent intermolecular modes. We note that the intramolecular modes are generally more harmonic than intermolecular modes and the total energy of the mode is expected to be more equally partitioned between the kinetic and potential energies in the intramolecular case. Therefore, the present analysis based only on KE could be improved by supplementary potential energy analysis for the intermolecular modes.7 Since the preceding analysis shows the importance of the HOD bending mode in the early IVR process in the subpicosecond time scale, the KE of this mode is likely to show a strong dependence on the initial H-bond number. This can be seen from Figure 10, where the band KE of the HOD mode was resolved by the initial number of H-bonds. It clearly shows a significantly faster and larger energy rise in the weakly H-bonded subensemble compared to strongly H-bonded ones. This confirms that the H-bond environment of the HOD molecule strongly influences the effectiveness of the IVR to the HOD bending mode.

Figure 8. Changes in kinetic energies of several important bands: (a) OH stretch (OH-s) bands of the HOD and solvent H2O molecules. (b) Bending bands of HOD, solvent H2O, and their sum. (c) Intermolecular bands below 1200 cm−1. The band kinetic energies were obtained by integrating the nonequilibrium fKESDs over the frequency ranges given in Table 2 and, if necessary, taking their differences. The bond excitation energy was 7.660 kcal/mol. The average equilibrium value was subtracted from each curve.

significant time-dependent changes. From Figure 8a, the OH stretch mode of the HOD molecule is seen to play a minor role in the IVR process but the OH vibrations of the solvent H2O molecules take up more than 0.4 kcal/mol of KE (11% of the OD bond KE at T = 48 fs, 3.654 kcal/mol). (Note that nonzero starting values of the band energies are due to thermal fluctuations, and they introduce small uncertainty in the exact amount of transferred energy.) This kind of uphill energy transfer is likely to be overestimated in this study due to the classical description of vibration employed. In contrast, the KE of the HOD bending band in Figure 8b shows a rapid initial increase up to 0.7 ps then slowly decays to zero. The transient peak of this curve corresponds to ∼0.4 kcal/mol of energy transferred, more than one-half of the total transient absorption by the IVR process shown in Figure 4. The HOH bending band from the solvent also plays a significant role in the VER process, with transient absorption of more than 0.5 kcal/mol (14% of the initial OD bond KE) at 4 ps. Finally, the intermolecular bands below 1200 cm−1 from the solvent in Figure 8c begin to absorb KE after ∼1.5 ps and eventually take up ∼2.5 kcal/mol

5. CONCLUSIONS In this paper, we have developed a theoretical method to calculate an instantaneous KESD from nonequilibrium MD simulations, based on the notion of instantaneous power spectra. The method enables simultaneous spectral and temporal analyses of molecular kinetic energies without formal limitation on their precision. Since the method does not make J



Article

nonequilibrium processes. This aspect will be investigated more thoroughly in the future.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail; [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (No.: 20110020033 and 2013013156) and the Samsung Science and Technology Foundation (SSTF-1086-BA1401-12).

Figure 10. Changes in the kinetic energy of the HOD bending band for subensembles classified by the initial number of H-bonds (NHB) formed by the HOD molecule. The bond excitation energy was 7.660 kcal/mol. The average equilibrium value was subtracted from each curve.

■

REFERENCES

(1) McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976. (2) Kubo, R.; Toda, M.; Hashitsume, N. Statistical Physics II, Nonequilibrium Statistical Mechanics, 2nd ed.; Springer: Berlin, 1991. (3) Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and Quantum Corrections from Molecular Dynamics for Liquid Water. J. Chem. Phys. 1983, 79, 2375−2389. (4) Lobaugh, J.; Voth, G. A. A Quantum Model for Water: Equilibrium and Dynamical Properties. J. Chem. Phys. 1997, 106, 2400−2410. (5) Burnham, C. J.; Li, J.; Xantheas, S. S.; Leslie, M. The Parametrization of a Thole-Type All-Atom Polarizable Water Model from First Principles and Its Application to the Study of Water Clusters (N = 2−21) and the Phonon Spectrum of Ice Ih. J. Chem. Phys. 1999, 110, 4566−4581. (6) Yagasaki, T.; Saito, S. A Novel Method for Analyzing Energy Relaxation in Condensed Phases Using Nonequilibrium Molecular Dynamics Simulations: Application to the Energy Relaxation of Intermolecular Motions in Liquid Water. J. Chem. Phys. 2011, 134, 184503. (7) Jeon, J.; Cho, M. Redistribution of Carbonyl Stretch Mode Energy in Isolated and Solvated N-Methylacetamide: Kinetic Energy Spectral Density Analyses. J. Chem. Phys. 2011, 135, 214504. (8) Page, C. H. Instantaneous Power Spectra. J. Appl. Phys. 1952, 23, 103−106. (9) Lampard, D. G. Generalization of the Wiener−Khintchine Theorem to Nonstationary Processes. J. Appl. Phys. 1954, 25, 802− 803. (10) Yagasaki, T.; Saito, S. Energy Relaxation of Intermolecular Motions in Supercooled Water and Ice: A Molecular Dynamics Study. J. Chem. Phys. 2011, 135, 244511. (11) Steinel, T.; Asbury, J. B.; Zheng, J.; Fayer, M. D. Watching Hydrogen Bonds Break: A Transient Absorption Study of Water. J. Phys. Chem. A 2004, 108, 10957−10964. (12) Rezus, Y. L. A.; Bakker, H. J. On the Orientational Relaxation of HDO in Liquid Water. J. Chem. Phys. 2005, 123, 114502. (13) Rezus, Y. L. A.; Bakker, H. J. Orientational Dynamics of Isotopically Diluted H2O and D2O. J. Chem. Phys. 2006, 125, 144512. (14) Wang, Z.; Pang, Y.; Dlott, D. D. Hydrogen-Bond Disruption by Vibrational Excitations in Water. J. Phys. Chem. A 2007, 111, 3196− 3208. (15) Bakker, H. J.; Skinner, J. L. Vibrational Spectroscopy as a Probe of Structure and Dynamics in Liquid Water. Chem. Rev. 2009, 110, 1498−1517. (16) Fujisaki, H.; Zhang, Y.; Straub, J. E. Non-Markovian Theory of Vibrational Energy Relaxation and Its Applications to Biomolecular Systems. Adv. Chem. Phys. 2011, 145, 1−33. (17) Seabra, G. d. M.; Walker, R. C.; Elstner, M.; Case, D. A.; Roitberg, A. E. Implementation of the SCC-DFTB Method for Hybrid QM/MM Simulations within the Amber Molecular Dynamics Package. J. Phys. Chem. A 2007, 111, 5655−5664.

explicit use of the potential energy information, it greatly simplifies the spectral analysis compared to conventional approaches based on the normal-mode analysis. The method is also inherently dynamic in nature, and all dynamic processes taking place in the system are reflected in the calculated instantaneous KESD. As an application of the method, we have studied the vibrational relaxation of an excited OD stretch mode in HOD dissolved in water. For quantitative comparison, we have also performed an IR pump−probe experiment on the same system. The OD bond excitation and relaxation processes were simulated with a sudden initial increase in the atomic velocities of the bond and subsequent NEMD simulation employing a QM/MM potential for HOD/H2O molecules. The KESDs show time-dependent changes of all spectral bands under 4170 cm−1, and it was possible to determine the intra- and intermolecular energy transfer pathways by analyzing the time evolution of band kinetic energies. From this, we found that ∼20% of the excess KE is transiently transferred to the HOD bending mode within 1 ps, and the low frequency modes of the solvent takes up more than 75% of the energy in the long time limit. The OD bond vibrational relaxation was found to proceed with decay times of 0.12 ps (∼10%) and 2.42 ps (∼90%). The major decay time is ∼30% higher than the experimentally determined value of 1.794 ps. We attribute this difference to the uncertainties in the computational force field model, the presence of the short-time component found from computation, which is not well-resolved in experiment, and the lack of quantum effects18 in our description of vibrational relaxation. In addition, comparison of the KESD time profiles of the OD bond at different frequencies confirmed that the decay rate of the excited state population is nearly frequencyindependent once the classical anharmonicity effect is accounted for. This is consistent with the assumption often used in the modeling of the experimental pump−probe signal. Analysis of the instantaneous KESDs in terms of the initial Hbond number of the chromophore revealed that strong Hbonding significantly hinders the short-time intramolecular vibrational energy redistribution process with subpicosecond time scale that mainly involves the HOD bending mode. Experimental verification of this behavior and others related to the IVR process would require better experimental time resolution in the subpicosecond time scale. Finally, the two forms of the iKESD, fKESD and pKESD, carry a deeper connotation regarding the directionality of time correlations in K



Article

(18) Stock, G. Classical Simulation of Quantum Energy Flow in Biomolecules. Phys. Rev. Lett. 2009, 102, 118301. (19) Ernst, R. R.; Bodenhausen, G.; Wokaun, A. Nuclear Magnetic Resonance in One and Two Dimensions; Oxford University Press: Oxford, 1987. (20) Nguyen, P. H.; Stock, G. Nonequilibrium Molecular-Dynamics Study of the Vibrational Energy Relaxation of Peptides in Water. J. Chem. Phys. 2003, 119, 11350−11358. (21) Wu, Y.; Tepper, H. L.; Voth, G. A. Flexible Simple Point-Charge Water Model with Improved Liquid-State Properties. J. Chem. Phys. 2006, 124, 024503. (22) Jeon, J.; Cho, M. An Accurate Classical Simulation of a TwoDimensional Vibrational Spectrum: OD Stretch Spectrum of a Hydrated HOD Molecule. J. Phys. Chem. B 2014, 118, 8148−8161. (23) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M. et al. Amber 11; University of California: San Francisco, 2010. (24) Boucher, W. Azara 2.8; Department of Biochemistry, University of Cambridge: Cambridge, U.K., 2010. (25) Kim, S.; Kim, H.; Choi, J.-H.; Cho, M. Ion Aggregation in High Salt Solutions: Ion Network Versus Ion Cluster. J. Chem. Phys. 2014, 141, 124510. (26) Rey, R.; Møller, K. B.; Hynes, J. T. Hydrogen Bond Dynamics in Water and Ultrafast Infrared Spectroscopy. J. Phys. Chem. A 2002, 106, 11993−11996. (27) Lawrence, C. P.; Skinner, J. L. Ultrafast Infrared Spectroscopy Probes Hydrogen-Bonding Dynamics in Liquid Water. Chem. Phys. Lett. 2003, 369, 472−477. (28) Wang, Z.; Pakoulev, A.; Pang, Y.; Dlott, D. D. Vibrational Substructure in the OH Stretching Transition of Water and HOD. J. Phys. Chem. A 2004, 108, 9054−9063. (29) Walrafen, G. E. Raman and Infrared Spectral Investigations of Water Structure. In The Physics and Physical Chemistry of Water; Franks, F., Ed.; Springer: New York, 1972; Vol. 1; pp 151−214.

L


Simultaneous Spectral and Temporal Analyses of Kinetic Energies in

Recommend Documents