A Mixed Quantum-Classical Molecular Dynamics Study of the

Article pubs.acs.org/JPCB

A Mixed Quantum-Classical Molecular Dynamics Study of the Hydroxyl Stretch in Methanol/Carbon Tetrachloride Mixtures III: Nonequilibrium Hydrogen-Bond Dynamics and Infrared Pump− Probe Spectra Kijeong Kwac and Eitan Geva* Department of Chemistry, University of Michigan, Ann Arbor, Michigan 48109-1055, United States ABSTRACT: We present a mixed quantum-classical molecular dynamics study of the nonequilibrium hydrogen-bond dynamics following vibrational energy relaxation of the hydroxyl stretch in a 10 mol % methanol/carbon tetrachloride mixture and pure methanol. The ground and first-excited energy levels and wave functions are identified with the eigenvalues and eigenfunctions of the hydroxyl’s adiabatic Hamiltonian and as such depend parametrically on the configuration of the remaining, classically treated, degrees of freedom. The dynamics of the classical degrees of freedom are in turn governed by forces obtained by taking the expectation value of the force with respect to the ground or excited vibrational wave functions. Polarizable force fields and nonlinear mapping relations between the hydroxyl transition frequencies and dipole moments and the electric field along the hydroxyl bond are used, which were previously shown to quantitatively reproduce the experimental infrared steady-state absorption spectra and excited state lifetime [Kwac, K.; Geva, E. J. Phys. Chem. B 2011, 115, 9184; 2012, 116, 2856]. The relaxation from the firstexcited state to the ground state is treated as a nonadiabatic transition. Within the mixed quantum-classical treatment, relaxation from the excited state to the ground state is accompanied by a momentum-jump in the classical degrees of freedom, which is in turn dictated by the nonadiabatic coupling vector. We find that the momentum jump leads to breaking of hydrogen bonds involving the relaxing hydroxyl, thereby blue-shifting the transition frequency by more than the Stokes shift between the steadystate emission and absorption spectra. The subsequent nonequilibrium relaxation toward equilibrium on the ground state potential energy surface is thereby accompanied by red shifting of the transition frequency. The signature of this nonequilibrium relaxation process on the pump−probe spectrum is analyzed in detail. The calculated pump−probe spectrum is found to be in reasonable agreement with experiment, thereby providing further credibility to the underlying force fields and mixed quantumclassical methodology.

I. INTRODUCTION Hydrogen bonding (H-bonding) is one of the most important interactions in chemical and biological systems. Alcohols can form H-bonds due to their hydroxyl headgroup. At the same time, alcohols are also soluble in apolar solvents due to their aliphatic tail group. The H-bonding structure and dynamics of methanol, which is the simplest alcohol, have been studied extensively experimentally, theoretically, and computationally.1−36 Methanol/CCl4 mixtures represent a popular model system for performing such studies. The wealth of experimental and computational studies suggests that methanol form Hbonded oligomers in such mixtures and that the structure and dynamics of these oligomers vary with concentration. The gas-phase hydroxyl stretch frequency is 3681.5 cm−1 for the OH stretch of a methanol (MeOH) molecule and 2717.6 cm−1 for the OD stretch of a methanol-d molecule (MeOD).37 Upon forming H-bonds, the hydroxyl stretch band red-shifts and broadens, thereby making it a sensitive probe of H-bond structure and dynamics. Infrared (IR) spectroscopy of the © 2013 American Chemical Society

hydroxyl stretch is therefore a sensitive and powerful probe of H-bonding.15,16,21−31 The H-bond structure and dynamics in methanol/CCl4 mixtures, in particular, have been studied extensively via one-dimensional (1D) and two-dimensional (2D) IR spectroscopy of the methanol’s hydroxyl stretch.23−31 Methanol molecules in such mixtures can be classified as α, β, γ, δ, or ε according to their H-bonding status as shown in Figure 1. More specifically, α corresponds to non-H-bonded methanol molecules, β corresponds to methanol molecules that serve as H-bond acceptors, γ corresponds to methanol molecules that serve as H-bond donors, δ corresponds to methanol molecules that serve simultaneously as donor and acceptor, and ε corresponds to methanol molecules that serve simultaneously as donor for one H-bond and acceptor for two H-bonds. Received: April 15, 2013 Revised: May 22, 2013 Published: May 28, 2013 7737

dx.doi.org/10.1021/jp403726t | J. Phys. Chem. B 2013, 117, 7737−7749

The Journal of Physical Chemistry B

Article

relaxation from the excited to ground state. More specifically, the kinetic model assumes two time scales for H-bond breaking, 200 fs and 2.5 ps, and two time scales for H-bond reconstruction, 7 and 20 ps. The fast and slow time scales were assumed to correspond to direct and indirect H-bond breaking and reconstruction mechanisms, respectively. Here, direct and indirect correspond to breaking and forming an Hbond which is either associated or not associated with the original relaxing OD stretch, respectively.25,26 Our goal in this paper is to explore the nonequilibrium Hbond dynamics that follows vibrational energy relaxation within the framework of a molecular mixed quantum-classical dynamical model. Importantly, within the mixed quantumclassical treatment, population relaxation between quantum states is accompanied by energy release to the classical degrees of freedom in the form of a momentum jump. Thus, another goal of this paper is to elucidate the ability of the mixed quantum-classical treatment to account for the highly nontrivial nonequilibrium solvent dynamics that can follow irradiative relaxation in liquid solutions. The remainder of the paper is organized as follows. In Section II, we outline the mixed quantum-classical approach and simulation techniques. The results are presented and discussed in Section III. Summary and concluding remarks are provided in Section IV. A more detailed description of the mapping relations and optical response theory used for calculating spectra is provided in two appendices.

Figure 1. Classification of hydroxyls as α, β, γ, δ, and ε according to their level of participation in H-bonding.

In this paper we focus on deuterated methanol (MeOD). We consider the H-bond structure and dynamics in 10 mol % MeOD/CCl4 mixture and pure MeOD. The choice of 10 mol % MeOD/CCl 4 mixture is motivated by the recent experimental studies by Fayer and co-workers on this system.23−27,29 The choice of pure MeOD is meant to elucidate the transferability of the conclusions regarding H-bond structure and dynamics from the mixture to pure MeOD. In the case of the 10 mol % MeOD/CCl4 mixture, the α and β methanol molecules give rise to a well-resolved narrow spectral band centered at 2690 cm−1, with fwhm of 20 cm−1, while the γ, δ, and ε methanol molecules give rise to a single broad band that consists of overlapping sub-bands corresponding to the individual species. More specifically, the γ methanol molecules give rise to a red-shifted broader band centered at 2600 cm−1 with fwhm of 80 cm−1, and the δ and ε methanol molecules give rise to an even broader and more red-shifted band centered at 2490 cm−1 with a fwhm of 150 cm−1.13 In the case of pure MeOD, the α and β populations are negligibly small, and the absorption spectrum consists of a single broad band that corresponds to contributions from the γ, δ, and ε species.13 Extensive measurements of 1D and 2D IR spectra in MeOD/ CCl4 mixtures were reported by Fayer and co-workers.23−27,29 Frequency-resolved OD stretch excited-state lifetimes measurements revealed a correlation between shorter lifetimes and lower OD stretch frequencies (i.e., stronger H-bonding).25 Similar observations have also been made in pump−probe studies on methanol and ethanol by Woutersen et al.38 and Laenen et al.39 In two previous papers we have been able to reproduce these observations by combining a mixed quantumclassical methodology with polarizable force fields.13,14 The purpose of the present paper is to extend the methodology developed in refs 13 and 14 to simulating the nonequilibrium H-bond dynamics that follows nonradiative relaxation of the OD stretch from its first excited state to the ground state. Fayer and co-workers have studied this process by using frequency-resolved pump−probe and 2D IR spectroscopies.26,29 Of particular importance for us here is a negative spectral feature that was observed to emerge at delay times longer than the excited state lifetime. Fayer and co-workers attributed this negative feature to the above-mentioned nonequilibrium process, wherein H-bonds first break and then partially reform as the system reaches equilibrium on the ground state potential energy surface (PES). They were also able to fit their results to a kinetic model that postulated these processes. Based on their kinetic model, those authors have argued that H-bonds at the low frequency side of the δ band are more likely to break. This means that the stronger H-bonds are more likely to be broken by the energy released in the

II. METHODS A. Mixed Quantum-Classical Molecular Dynamics. The mixed quantum-classical method employed here is similar to that described in refs 13 and 14 and will therefore only be briefly outlined with emphasis on new elements. We consider a MeOD/CCl4 mixture or pure MeOD where one of the hydroxyls, referred to as the tagged hydroxyl, is treated quantum-mechanically, while the remaining untagged degrees of freedom (DOF), referred to as the bath, are treated classically. For a given configuration of the bath DOF, Q, we define the vibrational energy levels and corresponding stationary wave functions of the hydroxyl stretch as the eigenvalues and eigenfunctions of the adiabatic Hamiltonian: Had(q ̂, p ̂ ; Q)Ψn(q; Q) = En(Q)Ψn(q; Q)

(1)

where Had(q ̂, p ̂ ; Q) = Kq(p ̂) + V (q ̂, Q)

(2)

Here, q̂ and p̂ are the vibrational coordinate and momentum operators of the tagged hydroxyl, n = 0, 1, 2, ... correspond to the ground, first-excited, second-excited, and so forth, vibrational states of the tagged hydroxyl stretch, Kq(p̂) is the vibrational kinetic energy of the tagged hydroxyl, and V(q̂, Q) is the overall potential energy. Within the mixed quantumclassical molecular dynamics (MD) simulation, the dynamics of the bath DOF, Q, is governed by the classical equation of motion. For each instantaneous bath configuration, Q, Ĥ q(q̂, p̂; Q) is diagonalized to obtain the corresponding adiabatic e ig e n v a l u e s ( e n e r g y l e v el s ) an d e i g e n f u n c t i o n s , {En(Q),Ψn(q;Q)} of the tagged hydroxyl stretch. These adiabatic eigenvalues and eigenfunctions are then used to generate the forces exerted by the quantum DOF on the classical DOF. The reader is referred to refs 13 and 14 for a 7738



Article

than the fundamental frequency in light of the anharmonicity of the vibrational potential. As can be seen from Figure 2b, the bond length increases according to the following trend 12 > 11 > 01 > 00, which can also be traced back to anharmonicity and the enhanced H-bonding at the excited state.14 It should be noted that using the mapping relations between the electric field along the OD stretch and the OD bond length involves approximating the forces in the following manner:

more detailed discussion and the explicit expressions of these forces. Nonlinear mapping relations can be established between the electric field along the tagged OD bond and the corresponding bond-length and stretch frequency (see Figure 2).13,14 It should

⟨ψi|F(q ̂, Q )|ψi ⟩ ≈ F(⟨ψi|q|̂ ψi ⟩, Q )

i = 0 or 1

(3)

⟨ψi|F(q ̂, Q )|ψi ⟩ + ⟨ψj|F(q ̂, Q )|ψj⟩ 2 ⎛ ⟨ψi|q|̂ ψi ⟩ + ⟨ψj|q|̂ ψj⟩ ⎞ ≈ F⎜⎜ , Q ⎟⎟ 2 ⎝ ⎠ (i , j) = (0, 1) or (1, 2)

(4)

Once established, one may take advantage of these relations in order to avoid explicit on-the-fly diagonalization of the instantaneous adiabatic Hamiltonian, thereby reducing the computational cost of the mixed quantum-classical simulation to that of a fully classical MD simulation. This is particularly important when attempting to simulate nonequilibrium MD, as in this paper, which is significantly more computationally costly. All of the results in this paper were obtained using the abovementioned mapping relations rather than via direct on-the-fly diagonalization. The simulations were carried out with the AMBER 10 program package40 whose source code was modified in order to accommodate features which are unique for the mixed quantum-classical MD. B. Population Relaxation. In a previous paper, we calculated the rate constant, k01, for the nonadiabatic transition from the first-excited to the ground state of the OD stretch.14 The calculation was carried out within the framework of Fermi’s golden rule. A description of population relaxation via a rate constant implies exponential rate kinetics. More specifically, starting with the OD stretch in its excited state at time t = 0, the probability that it will make the transition to the ground state at a later time t ≥ 0 is given by

Figure 2. (a) Correlation plots of the OD stretch frequency (obtained via on-the-fly diagonalization) vs the electric field along the OD bond; (b) correlation plots of the OD bond length (obtained from the expectation value of the OD stretch displacement) vs the electric field along the OD bond. The solid lines correspond to the mapping relations obtained by polynomial fits (see Appendix A). The four correlation plots and mapping relations correspond to the different PES's that govern the bath dynamics (see text for details).

be noted that four different mapping relations are shown in Figure 2 (denoted 00, 01, 11, and 12). These correspond to the different types of bath dynamics involved in calculating the pump−probe spectra within the mixed quantum-classical methodology (see Section II C). More specifically, 00 and 11 correspond to bath dynamics on the adiabatic PES's that correspond to the ground (00) and first-excited (11) vibrational states of the tagged OD stretch, while 01 and 12 correspond to bath dynamics on PES's given by the arithmetic averages of the ground and first-excited PES's (01) and firstand second-excited PES's (12). It should also be noted that in the case of 00, 11, and 01 the OD stretch frequency corresponds to the fundamental transition between the ground and first-excited state, ω10 = (E1 − E0)/ℏ, while in the case of 12 the stretch frequency corresponds to the overtone transition frequency between the first- and second-excited states, ω21 = (E2 − E1)/ℏ. Finally, the bond lengths in the 00 and 11 cases are given by ⟨ψ0|q̂|ψ0⟩ and ⟨ψ1|q̂|ψ1⟩, respectively, while in the 01 and 12 cases they are given by the corresponding arithmetic averages, that is, (⟨ψ0|q̂|ψ0⟩ + ⟨ψ1|q̂|ψ1⟩)/2, and (⟨ψ1|q̂|ψ1⟩ + ⟨ψ2|q̂|ψ2⟩)/2, respectively. As can be seen from Figure 2a, the 00, 11, and 01 mapping relations between the fundamental frequency and the electric field essentially coincide, although the likelihood for configurations with large electric fields increases from the ground to the excited state. This trend can be traced back to the fact that H-bonding is stronger in the excited state.14 The mapping relation for the 12 PES is shifted to lower frequency, which can be traced back to the fact that the overtone frequency is lower

P01(t ) = k10e−k10t

(5)

Thus, the probability for a transition from the excited to the ground state within the simulation time step δ is given by δP01(t ) =

∫t

t+δ

dt ′P01(t ′) = e−k10t (1 − e−k10δ)

(6)

Due to the separation of time scales between the high frequency OD stretch and the slow accepting bath modes, one can assume that the transition is instantaneous on the bath time scale, that is, that the bath configuration, Q, after the transition is the same as it was before the transition. Energy conservation therefore dictates that the transition is accompanied by a momentum jump in the bath DOF. Within the mixed quantumclassical framework, this momentum jump is given by:41 ΔP(Q) = M1/2{sgn(d10 ̅̂ (Q) ·P ̅ ) ×

(d10 ̅̂ (Q) ̅̂ (Q) ·P ̅ )2 + 2ℏω01(Q) − (d10 ̅̂ (Q) ·P ̅ )}d10 (7) 1/2

1/2

Here, M is a diagonal matrix with (Mj) along the diagonal, P̅j = Pj/(Mj)1/2, d̅10,j = d10,j/(Mj)1/2, and d̅̂ = d̅10/|d̅10| is a unit 7739



Article

lifetime corresponds to the time scale of equilibration on the ground PES. Finally, we note that the pump−probe signal is usually reported as the real part of the Fourier−Laplace transform of Ppp(t2,t3) with respect to t3:

vector along d̅10. The vector d10 is the nonadiabatic coupling coefficient defined as d10 = ⟨ψ1(Q)|∇Q|ψ0(Q)⟩. It should be noted that one needs to explicitly diagonalize the adiabatic Hamiltonian at the time of the nonadiabatic transition in order to obtain the nonadiabatic coupling vector which is necessary for evaluating the momentum jump. However, since diagonalization is only required at a single point in time, rather than onthe-fly, this does not add significantly to the computational cost. C. Pump−Probe IR Spectra. The pump−probe signal can be expressed in terms of the third-order optical response function. The latter involves three field-matter interactions, separated by time intervals t1 (between the first and second interactions), t2 (between the second and third interactions) and t3 (between the third interaction and detection). In the case of impulsive pump−probe spectroscopy, t1 = 0, t2 is the time delay between the impulsive pump and probe pulses, and t3 is the time delay between the probe pulse and detection. Within the mixed quantum-classical Liouville formalism of optical response, one can show that the pump−probe signal is given by (see Appendix B): Ppp(t 2 , t3) =

⟨|μ010 |2 μt10 μt10+ t exp[− i 2 2 3

+ ⟨|μ010 |2 μt10 μt10+ t exp[− i 2

2

3

− ⟨|μ010 |2 μt21μt21+ t exp[− i 2

2

3

− ⟨|μ010 |2 μt10 μt10+ t exp[− i 2

2

3

∫t ∫t ∫t

t 2 + t3

∫t

t 2 + t3

Ppp(t 2 , ω3) = Re

2

dτω10(τ )]⟩00,00,01

(9)

(10)

We use QM/MM (based on the ONIOM method42 as implemented in the Gaussian 09 package43 and including all the molecules in the simulation box) in order to calculate μ′ for 100 randomly chosen configurations taken from equilibrium ground state mixed quantum-classical MD simulations. We establish a linear map that correlates μ′ with the electric field along the hydroxyl stretch. It should be noted that we have recently shown that accounting for the variation on μ′ considerably improves the agreement between the calculated absorption spectrum and the experimental one, particularly in the region of the high frequency αβ band.14 D. Simulation Details. We carried out mixed quantumclassical simulations at 300 K for the following two systems: (1) 10 mol % MeOD/CCl4 mixture consisting of 26 MeOD molecules and 230 CCl4 molecules in a cubical simulation cell, (2) pure methanol consisting of 256 MeOD molecules in a cubical simulation cell. Polarizable force fields were used as in ref 14. Periodic boundary conditions were used, and the simulation cell size was determined by the final value of a 800 ps classical trajectory under constant pressure (1.0 atm) and temperature (300 K) conditions. We start out by equilibrating over 200 ps using fully classical MD (via Langevin dynamics with collision frequency set to 1.0 ps−1 at 300 K). This was followed by sampling initial configurations which are 2 ps apart. For each chosen configuration, randomly chosen methanol molecules were tagged. The initial velocities were then randomized by reassigning them from a Maxwellian distribution at 300 K, followed by a 200 ps 300 K constant temperature fully classical MD simulation, using Langevin dynamics with the collision frequency of 1.0 ps−1. In the next step, a 35 ps 300 K constant temperature mixed quantum-classical trajectory on the 00 or 11 PES was simulated using the weak-coupling algorithm of Berendsen et al.44 Following equilibration, 50 ps mixed quantum-classical constant energy runs were carried out on the 00 or 11 PES’s. During these runs, configurations were sampled 2 ps apart that were used as initial configurations for nonequilibrium classical mapping MD production runs. Results obtained from two different types of nonequilibrium MD simulations are reported below:

2

dτω10(τ )]⟩00,11,01 dτω21(τ )]⟩00,11,12 dτω10(τ )]⟩00,(11,00),01

dt3eiω3t3{Ppp(t 2 , t3)}

μt jk = μ′⟨Ψj(Q t)|q|̂ Ψk(Q t)⟩

2

t 2 + t3

∞

The signal is therefore reported as a function of the delay time, t2, and detection frequency, ω3, where the latter is the Fourier frequency variable conjugated with the detection time variable t3. In other words, the pump−probe signal is reported as a t2dependent spectrum given as a function of ω3. It should be noted that the transition dipole moments in eq 8, {μjkt }, are treated as time-dependent due to their explicit dependence on the instantaneous bath configuration. In other words, non-Condon effects are accounted for. To this end, we follow the same procedure used in ref 14. More specifically, we assume that the transition dipole moment is given by

2

t 2 + t3

∫0

(8)

The four terms in eq 8 represent contributions from different Liouville pathways associated with distinctly different (t2 + t3)long nonequilibrium MD trajectories. All pathways start with equilibrium configurations sampled on the ground state, 00, PES. However, they differ with respect to the PES's used for the dynamics during the time period t2 + t3, as described below: • ⟨···⟩00,00,01 corresponds to propagation on the 00 PES during the time interval (0, t2), followed by propagation on the 01 PES during the time interval (t2, t2 + t3). • ⟨···⟩00,11,01 corresponds to propagation on the 11 PES during the time interval (0, t2), followed by propagation on the 01 PES during the time interval (t2, t2 + t3). • ⟨···⟩00,11,12 corresponds to propagation on the 11 PES during the time interval (0, t2), followed by propagation on the 12 PES during the time interval (t2, t2 + t3). • ⟨···⟩00,(11,00),01 corresponds to trajectories that start at time 0 on the 11 PES, but hop to the 00 PES within the time interval (0, t2) with probability P01(t2) (see eq 5). This is then followed by propagation on the 01 PES during the time interval (t2, t2 + t3). It should be noted that ⟨···⟩00,(11,00),01 → 0 when t2 ≪ k−1 10 , while ⟨···⟩00,11,01, ⟨···⟩00,11,12 → 0 when ≫ k−1 10 . It should also be noted that although the ⟨···⟩00,(11,00),01 term has a negative sign, it does not necessarily cancel out with the ⟨···⟩00,00,01 term. This is because the average over the nonequilibrium ensemble of bath configurations right after a vibrational energy relaxation event need not coincide with the ground state equilibrium ensemble of bath configurations. Thus, the lifetime of the signal is expected to be longer than the excited state lifetime, k−1 01 , and the time difference between the excited-state lifetime and signal 7740



Article

(1) The first type is designed to examine the effects of the momentum jump on the nonequilibrium dynamics following vibrational energy relaxation. In this case, we start with equilibrium configurations on the 11 PES. We then assume that starting at t = 1.0 ps, a nonadiabatic transition to the 00 PES can occur with probability P01(t − t0), where t0 = 1.0 ps. Importantly, a vibrational energy relaxation event is accompanied by a momentum jump, which can lead to structural dynamics in the vicinity of the relaxing OD stretch, and is followed by equilibration on the 00 surface until the system reaches thermal equilibrium on this PES. (2) The second type is designed to produce the pump− probe signal. In this case all trajectories start with equilibrium configurations on the 00 PES, followed by the four different (t2 + t3)-long nonequilibrium trajectories, as specified in eq 8. Importantly, the relaxation during t2 from the 11 PES to the 00 PES in the case of 00,(11,00),01 trajectory is treated as a stochastic process where the value of k01 depends on the Hbonding status of the tagged hydroxyl. More specifically, the values of k−1 10 used in the simulation is 2.15 ps, 2.15 ps, 1.0 ps, 0.5 ps, and 0.5 ps for α and β,23 γ,29 δ, and ε27 species, respectively. The time step used in the simulations is 1.0 fs. We considered 51 different values of t2 from 0 to 30 ps, equally spaced on a log scale. The following criteria were used to define an H-bond: (1) O···D distance is less than 2.65 Ǻ and (2) O···D−O angle is greater than 120° (”···” denotes the H-bond).

Figure 3. Average distance from the hydroxyl hydrogen of the tagged methanol vs the average magnitude of momentum change due to a momentum jump.

deuterium and untagged accepting oxygen is larger than that deposited into the tagged donating oxygen (note that no momentum is deposited into the untagged deuterium of a neighboring methanol) leads one to expect that H-bond breaking is more likely to turn a δ MeOD into a β MeOD, rather than into a γ MeOD (whereas H-bond breaking at the hydrogen site of a δ MeOD gives a β MeOD, H-bond breaking at the oxygen site gives a γ MeOD). B. Time Dependence of OD Subpopulations. Figure 4 shows the time dependence of the populations of the different

III. RESULTS A. Momentum Jump Acceptors. We start with analysis of nonequilibrium momentum jump simulations. In these simulations, each trajectory starts with the system in thermal equilibrium on the 11 PES’s. A vibrational energy relaxation event is then assumed to start at t = 1.0 ps. Emphasis is put on elucidating the nonequilibrium H-bond dynamics that follows vibrational energy relaxation. This dynamics occurs in two steps. The first step corresponds to breaking of H-bonds due to energy imparted by the momentum jump that accompanies the vibrational energy relaxation event. The second step corresponds to H-bond reconstruction as the system approaches its new equilibrium state on the 00 PES. The results reported in this section are based on averaging over 10 000 nonequilibrium trajectories in the case of the 10 mol % MeOD/CCl4 mixture and pure MeOD. Equation 7 dictates how the momentum jump is distributed among the bath DOF. To gain a better understanding of this distribution, we considered the dependence of the amount of momentum change of each bath atom on its distance from the deuterium of the relaxing OD stretch. The results for the 10 mol % MeOD/CCl4 mixture are shown in Figure 3 (the results for pure MeOD are very similar and therefore not shown). As can be seen, the momentum jump is concentrated on four atoms, namely, the deuterium of the relaxing OD stretch (at distance 0 Å), the oxygen of the relaxing OD stretch (at distance ∼1.0 Å), and the oxygen of the neighboring MeOD molecule that serves as H-bond acceptor in the case of Hbonding (at distance ∼1.8 Å) or the methyl carbon of the relaxing MeOD in the absence of H-bonding. Thus, when the relaxing OD is H-bonded, the momentum jump is localized on the deuterium and the donor and acceptor oxygens. This implies that H-bond breaking is most likely to involve an Hbond that the relaxing OD participates in. Furthermore, the fact that the amount of momentum deposited into the tagged

Figure 4. Time dependence of the populations of the different types of OD stretches (α, β, γ, δ) as the system undergoes vibrational energy relaxation (the small ε subpopulation is included in δ subpopulation). Vibrational energy relaxation starts at t = 1 ps. Parts a and b show the δ subpopulation while c and d show the α, β, and γ subpopulations. Parts a and c are for the 10 mol % MeOD/CCl4 mixture in CCl4, and b and d are for pure MeOD.

types of OD stretches (α, β, γ, δ) as the system undergoes vibrational energy relaxation (the ε subpopulation was added to the δ subpopulation). The populations at t < 1.0 ps fluctuate around their equilibrium values on the first-excited state (11) PES. Once vibrational energy relaxation starts at t = 1.0 ps, one observes a decrease of the δ subpopulation and an increase of the α, β, and γ subpopulations, within a time window of ∼2 ps after t = 1.0 ps. This trend suggests that the energy released as a momentum jump leads to breaking of H-bonds. The somewhat less pronounced increase in the γ subpopulation can be traced back to several factors. First, H-bond breaking can lead to both loss and gain of γ species. Second, as mentioned above, H-bond 7741



Article

observations reinforce the above-mentioned observation that breaking the H-bond on the oxygen site, δ → γ, is less ef ficient than breaking the H-bond on the hydrogen site, δ → β. We next consider the distribution of the delay times, trecomb − tbreak, between the momentum-jump-induced H-bond breaking time (tbreak) and the H-bond reforming time (trecomb) for the 10 mol % MeOD/CCl4 mixture and pure MeOD. The results are shown in Figure 6 for the case where the relaxing MeOD

breaking is more likely at the hydrogen site than at the oxygen site. The momentum-jump-induced H-bond breaking in the first step is followed by H-bond reconstruction as the system approaches thermal equilibrium on the ground state (00) PES during the second step. Equilibration is seen to occur on a somewhat slower time scale of ∼10 ps. It should be noted that H-bonding is only partially restored compared to its state prior to vibrational energy relaxation. This is because equilibrium on the 00 PES corresponds to less H-bonding in comparison to equilibrium on the 11 PES.14 Interestingly, the momentum-jump-induced relative drop in the δ population is larger for the 10 mol % mixture than for the pure methanol, even though the δ population in pure methanol is larger. More specifically, the δ population drops by 12.9% for the 10 mol % MeOD mixture and 6.4% for pure MeOD. The same trend is also reflected by the fact that the increase in α and β populations in the 10 mol % MeOD/CCl4 mixture are larger compared to these in pure MeOD. This difference is interpreted as evidence for the higher fragility of the incomplete H-bonding network in 10 mol % MeOD/CCl4 mixture in comparison to pure MeOD. C. Time Scales of H-Bond Breaking and Reforming. In this section we consider the time scales for the elementary Hbond breaking and reforming events. To this end, we considered the distribution of the delay times, tbreak − tjump, between the momentum jump time (tjump) and H-bond breaking time (tbreak) in the case where the tagged methanol belongs to the δ subpopulation at the time of the momentum jump. The results are shown in Figure 5 for 10 mol % MeOD/ CCl4 mixture (top panel) and pure MeOD (botttom panel).

Figure 6. Distribution of trecomb − tbreak for the same subset of molecules as in Figure 5. The molecules start as either β or γ (right after H-bond breaking) and turn into δ upon H-bond reconstruction. (a) 10 mol % MeOD/CCl4 mixture and (b) pure MeOD.

belongs to the δ subpopulation at the time of the momentum jump (same subpopulation as for Figure 5). The majority (77.8% and 84.8% for the 10 mol % MeOD/CCl4 and the pure MeOD, respectively) of the broken H-bonds are observed to undergo rapid recombination within ∼200 fs after being broken. Thus, while ∼30−40% δ H-bonds are broken by the momentum jump, more than two-thirds are restored within ∼200 fs, so that the net change of the population of δ species is ∼10% (see Figure 4). Finally, we examine the correlation between tbreak − tjump and trecomb − tbreak. A correlation plot of these two times is shown in Figure 7 for the case when the OD stretch belongs to the δ subpopulation at the time of the momentum jump (i.e., for the same subset of molecules as in Figures 5 and 6). For the δ → β → δ pathway (top panel in Figure 7), we observes an inverse correlation between the breaking and reconstruction times, so that the faster H-bonds to break are also the slower ones to reconstruct. This trend is attributed to the fact that an early Hbonding breaking event implies that the energy release associated with the momentum jump will take longer to dissipate, which will in turn slow down reconstruction. This should be contrasted with the δ → γ → δ pathway (bottom panel in Figure 7), where reconstruction is observed to be fast regardless of how fast the H-bond was broken. This trend can be attributed to the fact that the δ → γ pathway is less likely and therefore correspond to late H-bond breaking events. In this case, reconstruction does not require for the energy to first dissipate away. As we have seen above, the simulation results suggest that the majority of broken H-bonds reconstruct on a very short time scale (∼200 fs). Nevertheless, a smaller, but significant, subset of H-bonds take longer to reconstruct. We find that the

Figure 5. Distribution of tbreak − tjump for MeOD molecules that start out as δ and upon H-bond breaking turn into either β or γ. (a) 10 mol % MeOD/CCl4 mixture and (b) pure MeOD.

We first note that ∼30−40% of δ H-bonds break as a result of the momentum jump. Generally speaking, H-bond are observed to break within ∼200 fs following the momentum jump. There is little difference in H-bond breaking dynamics between the 10 mol % MeOD mixture and the pure MeOD. Breaking an H-bond implies that a δ MeOD turns into either a β MeOD or a γ MeOD. A δ MeOD is more likely to turn into a β MeOD (71.7%) than a γ MeOD (28.3%). Furthermore, producing γ is observed to take longer than producing β (as evidenced by the later emergence of the maxima). Both 7742



Article

Figure 8. (a) Nonequilibrium dynamics of the average OD stretch frequency during vibrational energy relaxation. Vibrational energy relaxation starts at t = 1.0 ps and is accompanied by a momentum jump. The latter leads to breaking of H-bonds, which is associated with the sharp initial blue-shift. This is followed by a red-shift, which is associated with partial H-bond reconstruction as the system reaches equilibrium on the ground-state PES. (b) Same as a, except that the plot for pure MeOD was shifted upward by 45 cm−1 in order to make it easier to compare the time dependence to that in the case of the 10 mol % MeOD/CCl4 mixture. Figure 7. Correlation between tbreak − tjump and trecomb − tbreak for the same subset of molecules as in Figures 5 and 6. (a) δ → β → δ, (b) δ → γ → δ.

t ∼ 2 ps. The overall blue shift associated with H-bond breaking following vibrational energy relaxation is ∼50 cm−1. This is somewhat larger than the Stokes shift (∼40 cm−1). The subsequent equilibration on the 00 PES is therefore accompanied by a red shift. This corresponds to the first scenario mentioned above and is consistent with experiment. Fitting the relaxation profile associated with the red shift in the second step to an exponential of the form e−t/τr yields 5.7 and 3.3 ps for τr in the cases of the 10 mol % MeOD/CCl4 mixture and the pure MeOD, respectively. The faster relaxation in pure MeOD can be attributed to the fact that having more methanol molecules around the tagged methanol accelerates the restoration of the broken H-bonds. The simulated time scales for the two subsequent steps associated with H-bond breaking and reconstruction can be compared to time scales for the same processes previously reported by Gaffney et al.25 and Asbury et al.26 The latter were obtained by fitting a kinetic model to the experimental pump− probe spectra. The kinetic model contained two rate constants for H-bond breaking, which were attributed to a fast and direct −1 pathway (k−1 bf = 200 fs) and a slow and indirect pathway (kbs = 2.5 ps) for H-bond breaking. The indirect H-bond pathway was conjectured to originate from breaking an H-bond associated with a methanol molecule other than the relaxing one and therefore requires energy transfer from the relaxing OD stretch to a neighboring OD stretch.25 The model also assumes two rate constants for H-bond reconstruction, which were −1 attributed to fast (k−1 f = 7 ps) and slow (ks = 20 ps) pathways. Our simulation results are consistent with the observation that H-bond breaking is significantly more rapid than H-bond reconstruction, although the experimental reconstruction time is somewhat longer. Our results are also consistent with the observation of two time scales of H-bond breaking and reconstruction, although not with the interpretation that they

majority of these follow the δ → β → δ pathway in the case of the 10 mol % MeOD/CCl4 mixture. D. Time Dependence of the Average OD Stretch Frequency. Breaking H-bonds as a result of the momentum jump is accompanied by a blue shift of the OD stretch frequency. If this blue shift is larger than the Stokes shift between the average steady state absorption and emission frequencies, it would be followed by a red shift as the OD stretch reaches equilibrium on the ground state PES. However, if this blue shift is smaller than the Stokes shift, it would be followed by an additional blue shift as the OD stretch reaches equilibrium on the ground state PES. Importantly, the first scenario would give rise to a distinctive negative feature in the pump−probe and two-dimensional IR spectra in the time window that corresponds to the equilibration on the groundstate PES. The fact that such a negative peak was observed experimentally26,29 lends credibility to the first scenario. We will now show that our molecular model and mixed quantumclassical treatment of the underlying dynamics are indeed able to capture this negative peak quantitatively. The time dependence of the average OD stretch frequency following momentum jump is shown in Figure 8 for the 10 mol % MeOD/CCl4 mixture and the pure MeOD. The average stretch frequency before t = 1 ps correspond to its excited state equilibrium value and is given by 2506 cm−1 and 2460 cm−1 for the 10 mol % MeOD/CCl4 mixture and the pure MeOD, respectively. The time dependence of the average OD stretch frequency in the 10 mol % MeOD/CCl4 mixture and pure MeOD are almost identical. In both cases, momentum jump, starting at t = 1 ps, is followed by a sharp blue shift that reaches its maximum value at 7743



Article

true at t > 200 fs and is due to the fact that the majority of weaker H-bonds already recombined within the initial 200 fs time window. F. Pump−Probe Spectra. Finally, the calculated pump− probe spectra for the 10 mol % MeOD/CCl4 mixture and pure MeOD are shown in Figure 10. The pump−probe spectra for

arise from direct and indirect pathways. More specifically, we observe that the two pathways correspond to δ → β (fast) and δ → γ (slow) pathways in the case of H-bond breaking and β → δ (slow) and γ → δ (fast) pathways in the case of H-bond reconstruction. E. Correlation between OD Stretch Frequency and HBond Breaking Rate. In this section we consider the correlation between the OD stretch frequency and the rate of H-bond breaking. To this end, we considered the population of MeOD molecules whose OD stretch frequency at t = tjump is below 2660 cm−1, which correspond to the range of the δ band, and divided them into the following three groups (see Figure 9a):

Figure 9. (a) Distribution of the OD stretch frequency at t = tjump for 10 mol % MeOD/CCl4 mixture. We divide the molecules into three groups labels L, M, and H. L correspond to MeOD molecules whose OD frequency is less than 2448 cm−1, M to molecules whose OD frequency is between 2448 and 2520 cm−1, and H to molecules whose OD frequency is between (2520−2660)cm−1. (b) Distribution of tbreak − tjump for the subsets L, M, and H.

Figure 10. Simulated pump−probe spectra for (a) 10 mol % MeOD/ CCl4 and (b) pure MeOD.

the 10 mol % MeOD/CCl4 mixture and pure MeOD are similar. Their interpretation is facilitated by considering the behavior at t2 < 1.0 ps separately from that at t2 > 1.0 ps, where ∼1.0 ps correspond to the first excited state lifetime. At t2 < 1.0 ps, the spectrum consists of two features: (1) A negative feature between ∼2200 and 2400 cm−1 that arises from the 00,11,12 Liouville pathway; (2) A positive feature between ∼2400 and 2600 cm−1 that arises from the 00,11,10 and 00,00,10 Liouville pathways. At t2 > 1.0 ps, the contributions from the 00,11,12 and 00,11,10 Liouville pathways diminish due to vibrational energy relaxation and a new contribution arises from the 00, (11,00),10 Liouville pathway. As a result, the pump−probe spectrum at t2 > 1.0 ps is seen to consist of two features: (1) A positive feature between ∼2400 and 2600 cm−1 that arises from the 00,00,10 Liouville pathway and (2) a negative feature which is blue-shifted relative to the positive feature and arises from the 00,(11, 00),10 Liouville pathway. The emergence of the negative peak is associated with momentum-jump-induced Hbond breaking that accompanies vibrational energy relaxation. Its disappearance is due to a dynamical red-shift associated with ground-state equilibration, which results in overlap between the negative and the positive features, until they cancel each other out at ∼30 ps. The shallow negative high-frequency feature at t2 > 1.0 ps seen in Figure 10 is the result of the two phase spectral dynamics described in Figure 8. Importantly, this feature is

• Group L with ωOD < 2448 cm−1. • Group M with 2448 cm−1 < ωOD < 2520 cm−1. • Group H with 2520 cm−1 < ωOD < 2660 cm−1. We then calculated the distribution of the H-bond breaking time, tbreak − tjump, separately for groups L, M, and H (see Figure 9b). The results show that the majority of H-bonds involving high frequency OD stretches break faster than low frequency ones. Since a higher OD stretch frequency correlates with a weaker H-bond, this implies that weaker H-bonds break faster. However, a closer examination shows that the fraction of H-bonds that break within the 200 fs window is also sensitive to which group they belong to. More specifically, 63.4%, 74.1%, and 81.1% of the H-bonds break within 200 fs for groups L, M, and H, respectively. Thus, the majority of H-bonds that break on longer time scales are associated with low frequency and therefore strongly H-bonded OD stretches. This observation may be compared with one of the assumptions made in the kinetic model used by Fayer and co-workers to analyze the pump−probe and two-dimensional spectra of the 10 mol % MeOD/CCl4 mixture,26,29 according to which stronger Hbonds associated with low-frequency OD stretches, are more likely to break.26 Our results seem to indicate that this is only 7744



Article

consistent with a similar negative feature observed experimentally in the 10 mol % MeOD/CCl4 mixture.26 It should be noted that this feature represents a signature of the nonequilibrium nature of the dynamics underlying pump−probe experiments, and therefore could not be reproduced via the standard methodology of computing nonlinear IR spectra from equilibrium ground state MD simulations. To the best of our knowledge this is the first time that such nonequilibrium effects were captured by a computational approach for modeling nonlinear IR spectra. Furthermore, the fact that our molecular model and simulation methodology are able to quantitatively reproduce this unique nonequilibrium feature, in addition to many other properties of this system, suggests that it provides an accurate description for the complex underlying dynamics of this system. Finally, it should be noted that while H-bond breaking leads to an increase of the β subpopulation, which is manifested by the emergence of a high frequency peak at ∼2690 cm−1, this does not give rise to the emergence of a pronounced blueshifted peak in the pump−probe spectrum. The lack of such a peak can be traced back to the suppression of high-frequency spectral features by non-Condon effects.14

■

MeOD/CCl4 mixture is consistent with experiment. Importantly, the shallow negative high-frequency feature at t2 > 1.0 ps was reproduced. This feature represents a signature of the nonequilibrium nature of the dynamics underlying pump−probe experiments and therefore could not be reproduced via the standard methodology of computing nonlinear IR spectra from equilibrium ground state MD simulations. To the best of our knowledge, this is the first time that such experimentally measurable nonequilibrium effects were captured by a computational approach for modeling nonlinear IR spectra.

APPENDIX A: MAPPING RELATIONS BETWEEN THE TRANSITION FREQUENCY AND BOND-LENGTH AND ELECTRIC FIELD ALONG THE OD BOND The mapping relations between the electric field along the OD bond and the transition frequency and OD bond length are given below for the 01 and 12 PES's. The corresponding mapping relations in the 00 and 11 cases were reported in refs 13 and 14. The frequency ω is in cm−1, the bond length is in Å, and the electric field x is given in units of e/a20.

IV. CONCLUDING REMARKS In this paper, we reported a comprehensive analysis of the nonequilibrium H-bond dynamics that follows vibrational energy relaxation in MeOD/CCl4 mixtures and pure MeOD, as well as its spectroscopic signature, within the framework of a molecular mixed quantum-classical dynamical model. Our main conclusions can be summarized as follows: • The momentum jump associated with vibrational energy relaxation events within our mixed quantum-classical scheme leads to the breaking of H-bonds associated with the relaxing MeOD molecule. • Breaking the H-bond on the oxygen site, δ → γ, is less efficient than breaking the H-bond on the hydrogen site, δ → β. • The majority of broken H-bonds reconstruct on a very short time scale (∼200 fs). Nevertheless, a smaller, but significant, subset of H-bonds take longer to reconstruct. We find that the majority of these follow the δ → β → δ pathway in the case of the 10 mol % MeOD/CCl4 mixture. • Our simulation results appear to be consistent with the observation that H-bond breaking is significantly more rapid than H-bond reconstruction, although the experimental reconstruction time is somewhat longer. Our results are also consistent with the observation of two time scales of H-bond breaking and reconstruction, although not with the interpretation that they arise from direct and indirect pathways. More specifically, we observe that the two pathways correspond to δ → β (fast) and δ → γ (slow) pathways in the case of H-bond breaking and β → δ (slow) and γ → δ (fast) pathways in the case of H-bond reconstruction. • The momentum-jump-induced H-bond breaking blue shifts the OD band by more than the Stokes shift between the steady state absorption and emission spectra. It is followed by a dynamical red shift associated with equilibration on the ground state surface. • The main features of the pump−probe signal as obtained from our simulations in the case of the 10 mol %

ω10(x) = 2707.3665 − 3564.1316x − 137014.19x 2 + 2656322x 3 − 27475254x 4

(A1)

ω21(x) = 2591.899997041 − 3753.044186592x − 140315.46484538x 2 + 8958691.306891x 3 − 514155706.9099x 4 + 15061345794.06x 5 − 233012045386.8x 6 + 1770082429456x 7 − 5105741849963x 8

(A2)

01 rOD (x) = 0.987400 + 0.379633x + 10.5689x 2

− 207.758x 3 + 2366.41x 4

(A3)

(12) rOD (x) = 1.01407 + 0.400531x + 6.75785x 2

+ 172.647x 3 − 7993.82x 4 + 140052x 5 − 780322x 6 (A4)

■

APPENDIX B: OPTICAL RESPONSE FOR A SINGLE-MODE THREE-STATE SYSTEM IN THE PRESENCE OF POPULATION RELAXATION In this appendix, we outline a mixed quantum-classical formulation of third-order optical response for a single-mode three-state system in the presence of population relaxation. Our formulation is inspired by mixed quantum-classical Liouville (MQCL) dynamics.45−61 However, we have incorporated several physically motivated simplifications in order to facilitate the implementation and tractability. The third-order response function is given by:62 S(3)(t1 , t 2 , t3) =

⎛ i ⎞3 ⎜ ⎟ Tr[μ G ̂ (t3)MG(t 2)MG(t1)Mρeq̂ ] ⎝ℏ⎠ (B1)

Here, G(t) is the field-free time evolution superoperator, M· = [μ̂,·] is the dipole moment superoperator and ̂ (R, P)⟨0(R)| ρeq̂ = |0(R)⟩ρeq,0 7745

(B2)



Article

momenta of the bath DOF, and ρklW(Z,t) = ⟨k(R)|ρ̂W(R,P;t)| l(R)⟩. The system dipole moment operator, μ̂ , is given by:

where ρ̂eq,0(R,P) is the bath density operator that correspond to equilibrium of the adiabatic ground-state PES. We start out by performing a partial Wigner transform49 over the bath DOF of the overall density operator, ρ̂(t): ̂ (R , P ; t ) = ρW

⎛ 1 ⎞ ⎜ ⎟ ⎝ 2π ℏ ⎠

N

μ ̂ = −μ10 (R)(|0(R)⟩⟨1(R)| + |1(R)⟩⟨0(R)|) − μ21(R)(|1(R)⟩⟨2(R)| + |2(R)⟩⟨1(R)|)

∫ dR′eiP·R′/ℏ⟨R − R′/2|ρ (̂ t )|R

+ R′/2⟩

(B3)

Within the adiabatic representation, one can describe the state of the overall system by a 3 × 3 matrix, where each of the elements is a bath phase space density: ⎛ ρ 00 (Z , t ) ρ 01(Z , t ) ρ 02 (Z , t )⎞ W W ⎜ W ⎟ ⎜ 10 11 12 ρ ̂(t ) ≈ ⎜ ρW (Z , t ) ρW (Z , t ) ρW (Z , t ) ⎟⎟ ⎜ 20 ⎟ 21 22 ⎝ ρW (Z , t ) ρW (Z , t ) ρW (Z , t )⎠

(B4)

Here, 0, 1, and 2 correspond to the ground, first-excited, and second-excited adiabatic system states, {|0(R)⟩, |1(R)⟩, and | 2(R)⟩}, respectively, Z = {R,P} stand for the coordinates and

21

11 12 ⎛ μ10 ρ 10 − μ10 ρ 01 ⎞ − μ10 ρW00 − μ21ρW02 − μ21ρW01 μ10 ρW μ10 ρW W W ⎜ ⎟ ⎜ 10 00 10 11 21 20 10 01 10 10 21 21 21 12 10 02 21 22 21 11 ⎟ ̂ ̂ ≈ μ ρW − μ ρW + μ ρW μ ρW − μ ρW + μ ρW − μ ρW μ ρW + μ ρW − μ ρW Mρ ̂ = μρ̂ ̂ − ρμ ⎜ ⎟ ⎜ 21 10 ⎟ 10 21 21 11 10 20 21 22 21 12 21 21 μ ρW − μ ρW − μ ρW μ ρW − μ ρW ⎝ μ ρW − μ ρW ⎠

In terms of trajectories, operating with M splits the trajectories with relative weights that depend on the transition dipole moments as indicated by eq B6. For example, if a trajectory is at a state {Z, 00, F} before operating with M, then operating with M splits it, with equal probabilities, into two trajectories, namely: {Z, 01, −μ10F} or {Z, 01, −μ10F}. In the absence of nonadiabatic coupling terms, the effect of field-free evolution on a trajectory is to propagate Z(t) classically on the αβ PES. Furthermore, in the case where α ≠ β, a phase factor exp(−i∫ t0ωαβZ(t)) is accumulated.61 We also assume that nonadiabatic transitions from the first-excited to the ground state need only be accounted for during t2. The rate constant for a nonadiabatic transition between the excited and ground state, k01, is described by Fermis golden rule and simulated as described in Section II B.14,59 The scheme described above gives rise to a mixed quantumclassical expression for the third-order optical response function that consists of contributions from 6 different types of nonequilibrium trajectories:

first-excited state (11) to the ground-state (00) during t2. The six contributions are explicitly given by: (3) S10,00,10 (t1 , t 2 , t3) = ⟨μ010 μt10 μt10+ t μt10+ t + t 1

[exp(−i + exp(i

(3) + S10,20,21 (t1 , t 2 , t3)}

+

∫0 ∫0

t1

1

dτω10(τ ) − i

2

∫t +t 1

t1

dτω10(τ ) − i

1

2

t1+ t 2 + t3

1

dτω10(τ ))

2

t1+ t 2 + t3

∫t +t

3

dτω10(τ ))]⟩10,00,10

2

(B8)

− CC (3) S10,(11,00),10 (t1, t 2 , t3) = −⟨μ010 μt10 μt10+ t μt10+ t 1

[exp(− i + exp(i

(3) (3) + S10,(11,00),10 (t1 , t 2 , t3) + S10,(11,11),10 (t1 , t 2 , t3)

+

(B6)

from trajectories that undergo population relaxation from the

⎛ i ⎞3 (3) S(3)(t1 , t 2 , t3) = ⎜ ⎟ {S10,00,10 (t1 , t 2 , t3) ⎝ℏ⎠ (3) S10,(11,11),21 (t1 , t 2 , t3)

(B5)

Here, we assumed that μ and μ are real (which is the case when the adiabatic wave functions are defined as real) and harmonic oscillator selection rules (expected to be valid when anharmonic effects are small). The dynamics of the phase space densities {ραβ W (Z,t)} is obtained by averaging over a large number of trajectories. Along each trajectory, we specify not just Z(t) but also the specific phase-space density that the trajectory contributes to [αβ](t) and a cumulative factor, F(t), that accounts for the effect of interaction with the field and phase accumulation during periods when α ≠ β (see below). Operating with the dipole moment superoperator, M, on ρ̂ has the effect of reshuffling the populations and coherences with the relative weights depending on the transition dipole moments: 10

∫0 ∫0

t1

dτω10(τ ) − i

1

∫t +t 1

t1

dτω10(τ ) − i

2

t1+ t 2 + t3

1

2 + t3

dτω10(τ ))

2

t1+ t 2 + t3

∫t +t

1

dτω10(τ ))]⟩10,(11,00),10

2

(B9)

− CC

(3) S10,20,10 (t1 , t 2 , t3)

(3) S10,(11,11),10 (t1, t 2 , t3) = ⟨μ010 μt10 μt10+ t μt10+ t 1

(B7)

[exp(− i

Here, α1β1, α2β2, and α3β3 correspond to the contributions from trajectories that occur on the α1β1 PES during the time interval (0, t1), then the α2β2 PES during the time interval (t1, t1 + t2), and finally the α3β3 PES during the time interval (t1 + t2, t1 + t2 + t3). α1β1,(11,00),α3β3 corresponds to contributions

+ exp(i − CC 7746

∫0 ∫0

t1

dτω10(τ ) − i

1

∫t +t 1

t1

dτω10(τ ) − i

2

t1+ t 2 + t3

1

2 + t3

dτω10(τ ))

2

t1+ t 2 + t3

∫t +t

1

dτω10(τ ))]⟩10,(11,11),10

2

(B10) dx.doi.org/10.1021/jp403726t | J. Phys. Chem. B 2013, 117, 7737−7749


Article

(3) S10,(11,11),21 (t1, t 2 , t3) = −⟨μ010 μt10 μt21+ t μt21+ t 1

∫0

[exp(− i + exp(i

∫0

t1

dτω10(τ ) − i

2

t1+ t 2 + t3

∫t +t 1

t1

1

∫t +t 1

(3) (ta , tb , tc) = ⟨μ010 μt10 μt10+ t μt10+ t + t exp( −i Pnr

2 + t3

a

dτω21(τ ))

−i

2

t1+ t 2 + t3

dτω10(τ ) − i

1

ta + tb + tc

∫t +t a

dτω21(τ ))]⟩10,(11,11),21 (B11)

−i

(3) S10,20,10 (t1 , t 2 , t3) = ⟨μ010 μt21μt21+ t μt10+ t + t exp( −i 1

∫0 −i

dτω10(τ ) − i

∫t

t1+ t 2

1

2

1

2

∫t +t 1

dτω20(τ )

−i

a

(B12)

∫0 −i

dτω10(τ ) − i

∫t

1

2

1

2

−i

∫t +t 1

dτω20(τ )

a

dτω21(τ ))⟩10,20,21 − CC

∫0

∞

dt 3

∫0

∞

dt 2

∫0

∞

dt1S(3)(t1 , t 2 , t3)E(t − t3)E

(t − t3 − t 2)E(t − t3 − t 2 − t1)

■ ■ ■

(B14)

⎪

λ = a, b, c

The rephasing and nonrephasing signals can then be associated with the components of the third-order polarization with the over phases ϕr = −ϕa + ϕb + ϕc and ϕnr = −ϕa + ϕb + ϕc, respectively:

−i

ta + tb + tc

∫t +t a

a

a

ta + tb + tc

∫t +t a

a

ta + tb + tc

a

b

a

b

c

a

ta + tb + tc

a

dτω10(τ )

∫0

ta

dτω10(τ )

b

a

b

c

∫0

ta

dτω10(τ )

dτω10(τ ))⟩10,(11,00),10

b

a

∫t +t

ta

dτω10(τ ))⟩10,(11,11),10

− ⟨μ010 μt10 μt21+ t μt21+ t + t exp(i −i

∫0

b

a

∫t +t

c

dτω10(τ ))⟩10,00,10

− ⟨μ010 μt10 μt10+ t μt10+ t + t exp(i −i

b

b

+ ⟨μ010 μt10 μt10+ t μt10+ t + t exp(i −i

a

b

b

a

b

c

∫0

ta

dτω10(τ )

dτω21(τ ))⟩10,(11,11),21

c

∫0

ta

dτω10(τ )

dτω10(τ ))⟩10,(11,00),10

b

a

b

c

∫0

ta

dτω10(τ )

dτω21(τ ))⟩10,(11,11),21

b

(B17)

(B18)

AUTHOR INFORMATION

REFERENCES

(1) Jorgensen, W. L. Quantum and Statistical Mechanical Studies of Liquids. 7. Structure and Properties of Liquid Methanol. J. Am. Chem. Soc. 1980, 102, 543−549. (2) Palinkas, G.; Hawlicka, E.; Heinzinger, K. A Molecular Dynamics Study of Liquid Methanol with a Flexible Three-Site Model. J. Phys. Chem. 1987, 91, 4334−4341. (3) Haughney, M.; Ferrario, M.; McDonald, I. R. MolecularDynamics Simulation of Liquid Methanol. J. Phys. Chem. 1987, 91, 4934−4940. (4) Matsumoto, M.; Gubbins, K. E. Hydrogen-Bonding in Liquid Methanol. J. Chem. Phys. 1990, 93, 1981−1994. (5) Meyer zum Büschenfelde, D.; Staib, A. Vibrational Spectroscopy and Molecular Dynamics of Solvated Methanol Tetramers and Pentamers. Chem. Phys. 1998, 236, 253−261. (6) Curtiss, L. A. Molecular Orbital Studies of Methanol Polymers Using a Minimal Basis Set. J. Chem. Phys. 1977, 67, 1144−1149. (7) Mo, O.; Yanez, M.; Elguero, J. Cooperative E_ects in the Cyclic Trimer of Methanol - An Ab-Initio Molecular Orbital Study. J. Mol. Struct.: THEOCHEM 1994, 120, 73−81. (8) Dixon, J. R.; George, W. O.; Hossain, M. F.; Lewis, R.; Price, J. M. Hydrogen-Bonded Forms of Methanol. IR Spectra and Ab-Initio Calculations. J. Chem. Soc., Faraday Trans. 1997, 93, 3611−3618. (9) Ohno, K.; Shimoaka, T.; Akai, N.; Katsumoto, Y. Relationship between the Broad OH Stretching Band of Methanol and HydrogenBonding Patterns in the Liquid Phase. J. Phys. Chem. A 2008, 112, 7342−7348. (10) Staib, A.; Borgis, D. A Quantum Multi-Mode Molecular Dynamics Approach to the Vibrational Spectroscopy of Solvated Hydrogen-Bonded Complexes. Chem. Phys. Lett. 1997, 271, 232−240.

(B15)

b

b

ACKNOWLEDGMENTS This project was supported by the National Science Foundation through grant CHE-1111495.

⎪

a

a

The authors declare no competing financial interest.

⎧∈[eiϕλe−iωt + e−iϕλeiωt ] t ≤ t ≤ t + δ λ λ Eλ(t ) = ⎨ ⎩ 0 otherwise

Pr(3)(ta , tb , tc) = ⟨μ010 μt10 μt10+ t μt10+ t + t exp(i

dτω10(τ )

Notes

where E(t) is the driving radiation field. We assume that the field consists of three consecutive impulsive pulses centered at ta < tb < tc (assumed to have square envelope for the sake of simplicity):

a

c

ta

(3) (3) Ppp (tb , tc) = Pr(3)(ta = 0, tb , tc) ≡ Pnr (ta = 0, tb , tc)

The third-order polarization is given by: P(3)(t ) =

dτω10(τ )

The impulsive pump−probe signal in the time domain corresponds to the case where the first and second pulses coincide (see eq 9)

(B13)

2

b

a

ta + tb + tc

∫t +t

ta

b

a

1

t1+ t 2 + t3

b

− ⟨μ010 μt10 μt21+ t μt21+ t + t exp( −i

3

∫0

dτω10(τ ))⟩10,(11,11),10

a

ta + tb + tc

∫t +t

c

b

a

2

1

a

− ⟨μ010 μt10 μt10+ t μt10+ t + t exp( −i

dτω10(τ ))⟩10,20,10 − CC

t1+ t 2

b

∫0

3

(3) (t1 , t 2 , t3) = −⟨μ010 μt21μt21+ t μt10+ t + t exp( −i S10,20,21 t1

b

dτω10(τ ))⟩10,00,10

a

ta + tb + tc

∫t +t a

1

t1+ t 2 + t3

a

b

a

t1

b

+ ⟨μ010 μt10 μt10+ t μt10+ t + t exp( −i

2

− CC

a

(B16) 7747



Article

(11) Staib, A. A Theoretical Study of Hydrogen Bond Dynamics of Methanol in Solution. J. Chem. Phys. 1998, 108, 4554−4562. (12) Veldhuizen, R.; de Leeuw, S. W. Molecular Dynamics Study of the Thermodynamic and Structural Properties of Methanol and Polarizable/Non-Polarizable Carbon Tetrachloride Mixtures. J. Chem. Phys. 1996, 105, 2828−2836. (13) Kwac, K.; Geva, E. A Mixed Quantum-Classical Molecular Dynamics Study of the Hydroxyl Stretch in Methanol/CarbonTetrachloride Mixtures: Equilibrium Hydrogen-Bond Structure and Dynamics at the Ground State and the Infrared Absorption Spectrum. J. Phys. Chem. B 2011, 115, 9184−9194. (14) Kwac, K.; Geva, E. Mixed Quantum-Classical Molecular Dynamics Study of the Hydroxyl Stretch in Methanol/CarbonTetrachloride Mixtures II: Excited State Hydrogen Bonding Structure and Dynamics, Infrared Emission Spectrum, and Excited State Lifetime. J. Phys. Chem. B 2012, 116, 2856−2866. (15) Errera, J.; Mollet, P. Intermolecular Forces and O-H Absorption Bands in Alcohols at 3μ. Nature 1936, 138, 882. (16) Liddel, U.; Becker, E. D. Infrared Spectroscopic Studies of Hydrogen Bonding in Methanol, Ethanol and tert-Butanol. Spectrochim. Acta, Part A 1957, 10, 70−84. (17) Bellamy, L. J.; Pace, R. J. Hydrogen Bonding by Alcohols and Phenols I. Nature of Hydrogen Bond in Alcohol Dimers and Polymers. Spectrochim. Acta 1966, 22, 525−533. (18) Bonner, O. D. A Comparison of Hydrogen- and DeuteriumBonding in Carbon Tetrachloride Solutions of Methanol. J. Chem. Thermodyn. 1970, 2, 577−581. (19) Graener, H.; Ye, T. Q.; Laubereau, A. Ultrafast Vibrational Predissociation of Hydrogen Bonds: Mode Selective Infrared Photochemistry in Liquids. J. Chem. Phys. 1989, 91, 1043−1046. (20) Laenen, R.; Rauscher, C. Transient Hole-Burning Spectroscopy of Associated Ethanol Molecules in the Infrared: Structural Dynamics and Evidence for Energy Migration. J. Chem. Phys. 1997, 106, 8974− 8980. (21) Bertie, J. E.; Zhang, S. L. Infrared Intensities of Liquids XXI: Integrated Absorption Intensities of CH3OH,CH3OD and CD3OD and Dipole Moment Derivatives of Methanol. J. Mol. Struct. 1997, 413−114, 333−363. (22) Kristiansson, O. Investigation of the OH Stretching Vibration Of CD3OH in CCl4. J. Mol. Struct. 1999, 477, 105−111. (23) Levinger, N. E.; Davis, P. H.; Fayer, M. D. Vibrational Relaxation of the Free Terminal Hydroxyl Stretch in Methanol Oligomers: Indirect Pathway to Hydrogen Bond Breaking. J. Chem. Phys. 2001, 115, 9352−9360. (24) Gaffney, K. J.; Piletic, I. R.; Fayer, M. D. Hydrogen Bond Breaking and Reformation in Alcohol Oligomers Following Vibrational Relaxation of a Non-Hydrogen-Bond Donating Hydroxyl Stretch. J. Phys. Chem. A 2002, 106, 9428−9435. (25) Gaffney, K. J.; Davis, P. H.; Piletic, I. R.; Levinger, N. E.; Fayer, M. D. Hydrogen Bond Dissociation and Reformation in Methanol Oligomers Following Hydroxyl Stretch Relaxation. J. Phys. Chem. A 2002, 106, 12012−12023. (26) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Goun, A.; Fayer, M. D. Hydrogen Bond Dynamics Probed with Ultrafast Infrared Heterodyne-Detected Multidimensional Vibrational Stimulated Echoes. Phys. Rev. Lett. 2003, 91, 23742. (27) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Fayer, M. D. Hydrogen Bond Breaking Probed with Multidimensional Stimulated Vibrational Echo Correlation Spectroscopy. J. Chem. Phys. 2003, 119, 12981−12997. (28) Asbury, J. B.; Steinel, T.; Stromberg, C.; Gaffney, K. J.; Piletic, I. R.; Goun, A.; Fayer, M. D. Ultrafast Heterodyne Detected Infrared Multidimensional Vibrational Stimulated Echo Studies of Hydrogen Bond Dynamics. Chem. Phys. Lett. 2003, 374, 362−371. (29) Asbury, J. B.; Steinel, T.; Fayer, M. D. Hydrogen Bond Networks: Structure and Evolution after Hydrogen Bond Breaking. J. Phys. Chem. B 2004, 108, 6544−6554. (30) Asbury, J. B.; Steinel, T.; Stromberg, C.; Corcelli, S. A.; Lawrence, C. P.; Skinner, J. L.; Fayer, M. D. Water Dynamics:

Vibrational Echo Correlation Spectroscopy and Comparison to Molecular Dynamics Simulations. J. Phys. Chem. A 2004, 108, 1107−1119. (31) Asbury, J. B.; Steinel, T.; Fayer, M. D. Vibrational Echo Correlation Spectroscopy Probes of Hydrogen Bond Dynamics in Water and Methanol. J. Lumin. 2004, 107, 271−286. (32) Gulmen, T. S.; Sibert, E. L., III. Fluctuating Energy Level Landau-Teller Theory: Application to the Vibrational Energy Relaxation of Liquid Methanol. J. Phys. Chem. A 2005, 109, 5777− 5780. (33) Gulmen, T. S.; Sibert, E. L., III. Vibrational Energy Relaxation of the OH(D) Stretch Fundamental of Methanol in Carbon Tetrachloride. J. Chem. Phys. 2005, 123, 204508. (34) Iwaki, L. K.; Dlott, D. D. Three-Dimensional Spectroscopy of Vibrational Energy Relaxation in Liquid Methanol. J. Phys. Chem. A 2000, 104, 9101−9112. (35) Iwaki, L. K.; Dlott, D. D. Ultrafast Vibrational Energy Redistribution within C−H and O−H Stretching Modes of Liquid Methanol. Chem. Phys. Lett. 2000, 321, 419−425. (36) Wang, Z.; Pakoulev, A.; Dlott, D. D. Watching Vibrational Energy Transfer in Liquids with Atomic Spatial Resolution. Science 2002, 296, 2201−2203. (37) Serrallach, A.; Meyer, R.; Günthard, H. H. Methanol and Deuterated Species: Infrared Data, Valence Force Field, Rotamers, and Conformation. J. Mol. Spectrosc. 1974, 52, 94−129. (38) Woutersen, S.; Emmerichs, U.; Bakker, H. J. A Femtosecond Midinfrared Pump-Probe Study of Hydrogen-Bonding in Ethanol. J. Chem. Phys. 1997, 107, 1483−1490. (39) Laenen, R.; Gale, G. M.; Lascoux, N. IR Spectroscopy of Hydrogen-Bonded Methanol: Vibrational and Structural Relaxation on the Femtosecond Time Scale. J. Phys. Chem. A 1999, 103, 10708− 10712. (40) Case, D.; Darden, T.; Cheatham, T., III; Simmerling, C.; Wang, J.; Duke, R.; Luo, R.; Crowley, M.; Walker, R.; Zhang, W. et al. AMBER 10; University of California: San Francisco, CA, 2008. (41) Hammes-Schiffer, S.; Tully, J. C. Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions. J. Chem. Phys. 1994, 101, 4657−4667. (42) Dapprich, S.; Komaromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A New ONIOM Implementation in Gaussian 98. 1. The Calculation of Energies, Gradients and Vibrational Frequencies and Electric Field Derivatives. J. Mol. Struct.: THEOCHEM 1999, 462, 1− 21. (43) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09; Gaussian Inc.: Wallingford, CT, 2009. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (45) Martens, C. C.; Fang, J. Semiclassical-Limit Molecular Dynamics on Multiple Electronic Surfaces. J. Chem. Phys. 1997, 106, 4918. (46) Donoso, A.; Martens, C. C. Simulation of Coherent Nonadiabatic Dynamics Using Classical Trajectories. J. Chem. Phys. 1998, 102, 4291. (47) Donoso, A.; Martens, C. C. Semiclassical Multistate Liouville Dynamics in the Adiabatic Representation. J. Chem. Phys. 2000, 112, 3980. (48) Donoso, A.; Kohen, D.; Martens, C. C. Simulation of Nonadiabatic Wave Packet Interferometry Using Classical Trajectories. J. Chem. Phys. 2000, 112, 7345. (49) Schütte, C. Konard-Zuse-Zentrum für informationstechnik Berlin, June 1999; p SC 99-10. (50) Santer, M.; Manthe, U.; Stock, G. Quantum-Classical Liouville Description of Multi-Dimensional Nonadiabatic Molecular Dynamics. J. Chem. Phys. 2001, 114, 2001. (51) Kapral, R.; Ciccotti, G. Mixed Quantum Classical Dynamics. J. Chem. Phys. 1999, 110, 8919−8929. 7748



Article

(52) Nielsen, S.; Kapral, R.; Ciccotti, G. Mixed Quantum Classical Surface Hopping Dynamics. J. Chem. Phys. 2000, 112, 6543. (53) Wan, C.; Schofield, J. Mixed Quantum-Classical Molecular Dynamics: Aspects of Multithreads Algorithms. J. Chem. Phys. 2000, 113, 7047. (54) Wan, C.; Schofield, J. Exact And Asymptotic Solutions of the Mixed Quantum-Classical Liouville Equation. J. Chem. Phys. 2000, 112, 4447. (55) Kapral, R. Progress in the Theory of Mixed-Quantum Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129. (56) Shi, Q.; Geva, E. A Derivation of the Mixed Quantum-Classical Liouville Equation from the Inuence Functional Formalism. J. Chem. Phys. 2004, 121, 3393−3404. (57) Hanna, G.; Kapral, R. Quantum-Classical Liouville Dynamics of Nonadiabatic Proton Transfer. J. Chem. Phys. 2005, 122, 244505. (58) Hanna, G.; Kapral, R. Non-Adiabatic Dynamics of Condensed Phase Rate Processes. Acc. Chem. Res. 2006, 39, 21. (59) Hanna, G.; Geva, E. Vibrational Energy Relaxation of a Hydrogen-Bonded Complex Dissolved in a Polar Liquid via the Mixed Quantum-Classical Liouville Method. J. Phys. Chem. B 2008, 112, 4048−4058. (60) Hanna, G.; Kapral, R. Quantum-Classical Liouville Dynamics of Proton and Deuteron Transfer Rates in a Solvated Hydrogen Bonded Complex. J. Chem. Phys. 2008, 128, 164520. (61) Hanna, G.; Geva, E. Multi-Dimensional Spectra via the Mixed Quantum-Classical Liouville Method: Signatures of Nonequilibrium Dynamics. J. Phys. Chem. B 2009, 113, 9278. (62) Mukamel, S. Principles of Nonlinear Optical Spectroscopy; Oxford: New York, 1995.

7749


A Mixed Quantum-Classical Molecular Dynamics Study of the

Recommend Documents