Mixed Quantum-Classical Liouville Approach for Calculating Proton

Article pubs.acs.org/JCTC

Mixed Quantum-Classical Liouville Approach for Calculating ProtonCoupled Electron-Transfer Rate Constants Farnaz Shakib and Gabriel Hanna* Department of Chemistry, University of Alberta, Edmonton, Alberta T6G 2G2, Canada ABSTRACT: In this work, we derive a general mixed quantum-classical formula for calculating thermal proton-coupled electron-transfer (PCET) rate constants, starting from the time integral of the quantum flux−flux correlation function. This formula allows for the direct simulation of PCET reaction dynamics via the mixed quantum-classical Liouville approach. Owing to the general nature of the derivation, this formula does not rely on any prior mechanistic assumptions and can be applied across a wide range of electronic and protonic coupling regimes. To test the validity of this formula, we applied it to a reduced model of a condensed-phase PCET reaction. Good agreement with the numerically exact rate constant is obtained, demonstrating the accuracy of our formalism. We believe that this approach constitutes a solid foundation for future investigations of the rates and mechanisms of a wide range of PCET reactions.

1. INTRODUCTION Coupling of electron transfer to proton transfer, or protoncoupled electron transfer (PCET), constitutes an important class of charge-transfer reactions in chemistry and biology.1−7 Because of its role in vital life processes such as photosynthesis, respiration, and DNA repair, PCET has been the subject of intense experimental4,8−14 and theoretical3,8,15−22 research over the past two decades. In condensed-phase PCET, the motions of the transferring electrons and protons, donor−acceptor atoms, and surrounding solvent atoms are coupled to varying degrees and can occur on different timescales. Thus, a general theoretical approach for studying the rates and mechanisms of PCET reactions should account for the quantum natures of the transferring electrons and protons and for a wide range of coupling strengths and time scales. Over the past two decades, the theoretical work in this field has largely focused on understanding the rates of PCET reactions. Much of this work was motivated by Marcus theory23,24 for electron transfer and by theories for nonadiabatic proton transfer.25−27 A host of golden-rule-based rateconstant expressions have been derived (after a series of approximations) and used to study PCET reactions across a range of coupling regimes, which depend on the nonadiabaticity of the reaction; that is, reactions could be adiabatic or nonadiabatic with respect to electron transfer and proton transfer.3,15,18,28,29 In the 1990s, Cukier extended Marcus theory for electron transfer by introducing a proton-transfer coordinate.15,17,28 The resulting rate-constant expressions contain overlaps of the proton vibrational wave functions, which are localized in asymmetric double-well potentials for the reactant and product states. Following this, Soudackov and Hammes-Schiffer proposed a two-dimensional generalization of this theory, in which two collective solvent coordinates corresponding to the electron and proton transfer are introduced.18,30 Within this theory, the solvent environment © 2016 American Chemical Society

has been described by both a dielectric continuum model and an explicit molecular solvent−protein environment. Later, this theory was extended to incorporate the effects of the proton donor−acceptor vibrational motion.31−34 Analytical expressions have been derived for both concerted and sequential PCET reactions; however, it can be difficult to know which expression to use without prior mechanistic knowledge of the reaction under study. To evaluate these expressions, the input parameters (e.g., reorganization energies, equilibrium freeenergy differences, and couplings) must be defined in terms of quantities that can be determined computationally via electronic structure calculations and molecular dynamics simulations and/or experimentally. Rate theories that involve the direct simulation of PCET reaction dynamics are desirable because they provide a general approach for simultaneously studying the rates, mechanisms, and factors that affect PCET reactions across a wide range of nonadiabatic coupling regimes. Unlike the previous theories, such theories do not rely on prior mechanistic knowledge of the reaction. Previously, fewest switches surface hopping (FSSH)35 and ring-polymer molecular dynamics (RPMD)36,37 have been used to simulate PCET reaction dynamics.21,38−41 Of particular note, Kretchmer and Miller have extended the RPMD method to calculate thermal rate constants through direct simulations of the PCET reaction dynamics.21 RPMD is a highly efficient approximate quantum dynamics approach, but it can suffer from ergodicity issues and the neglect of real-time quantum coherence (because it is an imaginary-time path-integral-based approach). The latter becomes an issue when the quantum coherences persist beyond thermal time βℏ. Breakdowns have been reported for gas-phase scattering processes42 and for electron transfer in the inverted Marcus regime,43 which is Received: April 10, 2016 Published: May 27, 2016 3020

DOI: 10.1021/acs.jctc.6b00362 J. Chem. Theory Comput. 2016, 12, 3020−3029

Article

Journal of Chemical Theory and Computation

̂ where neq A is the equilibrium density of the reactant species, jA = (i/ℏ)[Ĥ ,Â ] is the flux of reactant species operator Â , jB̂ = (i/ ℏ)[Ĥ ,B̂ ] is the flux of product species operator B̂ , and the flux− flux correlation function is expressed in terms of the Kubotransformed correlation function

particularly relevant for PCET (although the breakdown in the inverted Marcus regime can be overcome by using the kinetically constrained RPMD method44). Thus, mixed quantum-classical and semiclassical methods, which incorporate real-time quantum coherences into the dynamics, can provide viable ways for simulating thermal and photoinduced PCET reactions in regimes where RPMD may break down. With regard to mixed quantum-classical methods, FSSH has been used to simulate the nonadiabatic dynamics of the classical degrees of freedom (DOF) involved in PCET reactions on mixed electronic−protonic adiabatic potential energy surfaces (PESs).38,40 However, FSSH suffers from an inability to treat decoherence, which can lead to inaccurate results.45 For example, it was shown that FSSH yields the incorrect scaling of the electron-transfer rate constant with the diabatic coupling46,47 as expected from Marcus’s golden rule expression, and the incorrect long-time dynamics of the spin-boson model in the strongly nonadiabatic, biased, and strong system−bath coupling regimes.48 Several attempts have been made at correcting this decoherence issue, but owing to the lack of a formal derivation, they are ad hoc in nature.49−52 Recently, we have demonstrated that the mixed quantum-classical Liouville (MQCL) approach53−62 captures the exact state population dynamics of model PCET reactions more accurately than does FSSH.63,64 This was attributed to the ability of the MQCL approach to rigorously treat quantum coherence effects through mean-surface evolution. In this article, we derive an expression for calculating rate constants of thermal PCET reactions, starting from the time integral of the quantum flux−flux correlation function of the reactant and product species operators.65 To arrive at this expression, the full quantum equilibrium structure is simplified greatly by invoking a high-temperature approximation for the classical DOF and replacing the quantum dynamics of the correlation function with MQCL dynamics. This expression is general and therefore could be applied across all nonadiabatic coupling regimes and does not rely on any prior mechanistic assumptions. This article comprises the following sections: We derive an expression for the time-dependent rate coefficient (from which the rate constant may be extracted) in the Derivation of the Time-Dependent Rate Coefficient section. This expression is applied to a reduced PCET model, which we describe in the PCET Model section. The details of simulating this expression are presented in the Simulation Details section. In the Results and Discussion section, we present our rate constant result and compare it to the exact value obtained using the quasi-adiabatic path-integral method.66 Concluding remarks are made in the Conclusions section. Some details of the derivation are provided in Appendices 1 and 2.

⟨jÂ ; jB̂ (t )⟩ =

1 β

∫0

β

⎞ ̂⎛ i ̂ dλ e λH ⎜ [Ĥ , Â ]⎟e−λH jB̂ (t ) ⎝ℏ ⎠

(2)

where β = 1/kT. For the purposes of a simulation, it is convenient to consider the time-dependent rate coefficient defined as the finite time integral of the flux−flux correlation function kAB(t ) =

1 nAeq

∫0

t

dt ′⟨jÂ ; jB̂ (t ′)⟩ =

1 βnAeq

i ̂ [B(t ), Â ] ℏ (3)

In anticipation of a mixed quantum-classical description of the system, one then partitions the system into a subsystem of interest (which will be treated quantum mechanically) and a bath (which will be treated classically) and performs a partial Wigner transform over the bath variables to ultimately arrive at the following expression for the time-dependent rate coefficient: kAB(t ) =

2 ∑ ℏβnAeq αα ′

∫ dX Im[B Wαα′(X , t )W Aα′α(X , 0)]

(4)

where W Aα ′ α(X , 0) =

∑ ∫ dX′AWα α ′(X′)W α ′α α ′ α(X′, X , 0) 1 1

1 1

α1α1′

(5)

In the above equations, X ≡ (Q, P) denotes the vector of bath positions (Q) and momenta (P), subscript W indicates that a partial Wigner transform has been taken, and the superscripts denote a representation in the adiabatic basis, defined by the solutions of time-independent Schrödinger equation Ĥ W(Q)|α; Q⟩ = Eα(Q)|α; Q⟩. In eq 5, the initial values of the matrix elements of spectral density W(X′, X, 0) are given by69 W α1′α1α ′ α(X ′, X , 0) =

1 (2π ℏ)2ν ZQ

∫ dZ dZ′e−(i/ℏ)(P·Z+P′·Z′)

Z Z′ |Q ′ − ⟩|α1; Q ′⟩ 2 2 Z′ −βĤ Z × ⟨α1′; Q ′|⟨Q ′ + |e |Q − ⟩|α ; Q ⟩ 2 2 × ⟨α′; Q |⟨Q +

(6)

where ν is the position−space dimension, and ZQ is the partition function. Because ⟨Q + Z/2|Q′ − Z′/2⟩ = δ(Q + Z/2 − (Q′ − Z′/2)) is a delta function, Q + Z/2 = Q′ − Z′/2 and Z′ = 2(Q′ − Q − Z/2). Therefore, eq 6 may be simplified to give

2. DERIVATION OF THE TIME-DEPENDENT RATE COEFFICIENT 2.1. General Expression. For completeness, we begin with a brief summary of the quantum mechanical rate theory presented in refs 67 and 68 because it forms the basis of our derivation of a mixed quantum-classical expression for calculating PCET rate constants in the next subsection. The rate constant for a transformation A ⇌ B may be calculated from the time integral of the quantum mechanical flux−flux correlation function:65 ∞ 1 kAB = eq ∫ dt ⟨jÂ ; jB̂ (t )⟩ nA 0 (1)

W α1′α1α ′ α(X ′, X , 0) = ×

1 e−(2i/ ℏ)P ′·(Q ′− Q ) (2π ℏ)2ν ZQ

∫ dZ e−(i/ℏ)(P−P′)·Z⟨α′; Q |α1; Q ′⟩

× ⟨α1′; Q ′|⟨2Q ′ − Q −

Z −βĤ Z |e |Q − ⟩|α ; Q ⟩ 2 2

(7)

2.2. Mixed Quantum-Classical Expression for PCET. To compute the rates of PCET reactions using eq 4, we must specify the forms of reactant and product species operators Â 3021


Article

Journal of Chemical Theory and Computation and B̂ , respectively. The operator associated with the PCET reactant state may be defined as Â = |ϕe1⟩⟨ϕe1|h( −qp)

W Aα ′ α(X , 0) = ×

(8)

′

cijα1cklα1⟨ϕej|⟨ϕpi||ϕe1⟩⟨ϕe1|h( −qp)|ϕpk ⟩|ϕel⟩

∑ ∑

∫ dZ e−(i/ℏ)(P−P′)·Z

× e−(2i/ ℏ)P ′·(Q ′− Q )

where is the donor electronic diabatic state (or state 1) and h(−qp) is a Heaviside step function of protonic coordinate qp [i.e., if qp > 0, then h(qp) = 1, whereas if qp < 0, then h(qp) = 0]. Similarly, the operator associated with the PCET reactant state may be defined as

Z −β ∑i(Pi2 /2Mi) Z |e |Q − ⟩⟨α′; Q |α1; Q ′⟩ 2 2 ̂ Z × ⟨α1′; Q ′|e−β(h /2)(2Q ′− Q − Z /2)|α2 ; 2Q ′ − Q − ⟩ 2 Z −β(h ̂/2)(Q − Z /2) Z × ⟨α2 ; 2Q ′ − Q − |e |α3; Q − ⟩ 2 2 Z × ⟨α3; Q − |α ; Q ⟩ (14) 2

× ⟨2Q ′ − Q −

(9)

where |ϕ2e ⟩ is the acceptor electronic diabatic state (or state 2). To evaluate matrix element AαW1α′1 mentioned in eq 5, we expand adiabatic wave function |α; Q⟩ in an orthonormal set of twoparticle basis functions as |α ; Q ⟩ =

∑ ∫ dX ′ α1α1′α2α3

i , k j , l = 1,2

|ϕ1e ⟩

B̂ = |ϕe2⟩⟨ϕe2|h(qp)

1 (2π ℏ)2ν ZQ

To further simplify WαA′α(X, 0), we next consider the ⟨2Q′ − Q 2

− Z/2|e−β∑i(Pi /2Mi)|Q − Z/2⟩ term in eq 14. In Appendix 1, we

α (Q )|ϕpm⟩|ϕen⟩ ∑ cmn

(10)

m,n

show that it simplifies to where |ϕmp ⟩ is chosen to be the mth eigenfunction of the quantum harmonic oscillator for a proton and index n runs over the electronic diabatic states 1 (donor) and 2 (acceptor). Consequently, we find that α1α1′ AW (Q ′)

=

⟨α1; Q ′||ϕe1⟩⟨ϕe1|h( − qp)|α1′;

=

∑ ∑

Z −β ∑i(Pi2 /2Mi) Z |e |Q − ⟩ 2 2 2 2 Mi e−(Mi /2βℏ )(2Q i − 2Q i′) 2 2βπ ℏ

⟨2Q ′ − Q − =

∏ i

Q ′⟩

′

cijα1cklα1⟨ϕej|⟨ϕpi||ϕe1⟩⟨ϕe1|h( − qp)|ϕpk ⟩|ϕel⟩

Inserting this expression back into eq 14 yields

i , k j , l = 1,2

(11)

We may then write W Aα ′ α(X , 0) = ×

WαA′α(X,

α1α1′

×e

W Aα ′ α(X , 0) =

0) in eq 5 as

×

1 (2π ℏ)2ν ZQ

∑ ∫ dX′∑ ∑

(15)

′


∑ ∑

∫ dZ e−(i/ℏ)(P−P′)·Z

× e−(2i/ ℏ)P ′·(Q ′− Q )

′

i , k j , l = 1,2

∫ dZ e

× ⟨α1′; Q ′|⟨2Q ′ − Q −

∑ ∫ dQ ′ dP ′ α1α1′α2α3

i , k j , l = 1,2

cijα1cklα1⟨ϕej|⟨ϕpi||ϕe1⟩⟨ϕe1|h( − qp)|ϕpk⟩|ϕel⟩

−(2i/ ℏ)P ′·(Q ′− Q )

1 (2π ℏ)2ν ZQ

−(i/ ℏ)(P − P ′)·Z

×

⟨α′; Q |α1; Q ′⟩

Z −βĤ Z |e |Q − ⟩|α ; Q ⟩ 2 2

∏ i

Mi

2

2

2βπ ℏ

2

e−(2Mi / βℏ )(Q i − Q i′)

× e−(β /2)Eα2(2Q ′− Q − Z /2) e−(β /2)Eα3(Q − Z /2) (12)

× ⟨α′; Q |α1; Q ′⟩⟨α1′; Q ′|α2 ; 2Q ′ − Q −

To simplify this expression, we first evaluate matrix element ̂ ⟨2Q′ − Q − Z/2|e−βH|Q − Z/2⟩. In Appendix 1, we show that after partitioning the Hamiltonian of the total system, Ĥ , into the kinetic energy of the environment and the remainder, ĥ, this matrix element reduces to

× ⟨α2 ; 2Q ′ − Q −

Z ⟩ 2

Z Z Z |α3; Q − ⟩⟨α3; Q − |α ; Q ⟩ 2 2 2 (16)

Because we are often interested in studying PCET reactions

̂ Z −βĤ Z |e |Q − ⟩ = e−β(h /2)(2Q ′− Q − Z /2) 2 2 ̂ Z −β ∑i(Pi2 /2Mi) Z × ⟨2Q ′ − Q − |e |Q − ⟩e−β(h /2)(Q − Z /2) 2 2

at high temperatures, we now invoke a high-temperature

⟨2Q ′ − Q −

approximation for the bath coordinates; that is, we assume that

(13)

lim

β→0

Substituting this expression back into eq 12 and inserting two completeness relations involving |α2⟩ and |α3⟩, we obtain

2

2

Mi /2βπ ℏ2 e−(2Mi / βℏ )(Q i − Q i′) = δ(Q i − Q i′). Inserting

these delta functions into WαA′α(X, 0) and integrating leads to 3022


Article

Journal of Chemical Theory and Computation W Aα ′ α(X , 0) =

1 (2π ℏ)2ν ZQ ′

∑ ∑ ∑

×

α1α1′α2α3 i , k j , l = 1,2

×

by directly simulating the PCET reaction dynamics using the MQCL approach, whereby the transferring electron and proton are treated quantum mechanically and the bath is treated classically. At t = 0, the rate coefficient is equal to zero because t the imaginary part of the phase factor, e(i/ℏ)∫0[Eα(t′)−Eα′(t′)] dt′, appearing in the MQCL-evolved matrix element is equal to zero.60,61 At early times, the rate coefficient rises to its transition-state theory value. If there is a sufficiently large separation between the time scales of the PCET reaction and the transient microscopic dynamics, then the rate coefficient decays to a plateau from which the rate constant, kAB, can be extracted (see Figure 1). In the long time limit, kAB(t) decays to zero.


∫ dP′ ∫ dZ e−(i/ℏ)(P−P′)·Z e−(β/2)E

α2(Q − Z /2)

× e−(β /2)Eα3(Q − Z /2)⟨α′; Q |α1; Q ⟩⟨α1′; Q |α2 ; Q −

Z ⟩ 2

Z Z Z |α3; Q − ⟩⟨α3; Q − |α ; Q ⟩ (17) 2 2 2 Next, using the definition of the Dirac delta function, ∫ dP′ e(i/ ℏ)P ′·Z = 2π ℏδ(Z), we carry out the integration over P′ to obtain × ⟨α2 ; Q −

W Aα ′ α(X , 0) =

1 (2π ℏ)ν ZQ ′

∑ ∑ ∑

×

α1α1′α2α3 i , k j , l = 1,2

×


∫ dZδ(Z)e−(i/ℏ)P·Z e−(β/2)E

α2(Q − Z /2)

× ⟨α′; Q |α1; Q ⟩⟨α1′; Q |α2 ; Q −

e−(β /2)Eα3(Q − Z /2) Figure 1. Schematic illustrating the time dependence of the rate coefficient kAB(t). The dotted line indicates the value of rate constant kAB.

Z ⟩ 2

Z Z Z |α3; Q − ⟩⟨α3; Q − |α ; Q ⟩ (18) 2 2 2 Carrying out the integration over Z and performing a series of simplifications yields the final expression for WαA′α(X, 0) (see Appendix 2 for details). × ⟨α2 ; Q −

W Aα ′ α(X , 0) =

e−βEα(Q ) ∑ ciα1 ′ckα1 (2π ℏ)ν ZQ i , k

3. PCET MODEL To demonstrate our mixed quantum-classical PCET rate theory, we applied it to a reduced model of a condensedphase PCET reaction presented by Ananth and Miller.66 This model is more complex than the model studied previously by the authors63,64 because the solvent polarization coordinate is now coupled to a bath of harmonic oscillators. This model has the following Hamiltonian:

0

∫−∞ dqpϕpi(qp)ϕpk(qp) (19)

Finally, the matrix element of the product species operator, ′ Bαα W (X, t), in eq 4 may be evaluated in a similar way to

̂ = HW

AαW1α1′(X′),

yielding the final expression for the time-dependent rate coefficient: 2 kAB(t ) = ∑ ℏβnAeq αα ′ ×

∫0

∞

∫

2mp

+

∑ j

P j2 2Mj

+

ps2 2ms

+ Vp(qp) + Ve(qp , qs)

+ Vps(qp , qs) + Vb(Q )

⎡ dX Im⎢∑ ciα2(t )ckα2′(t ) ⎢⎣ ik

(22)

where {qp, pp, mp}, {qs, ps, ms}, and {Qj, Pj, Mj} denote the positions, momenta, and masses of the proton, solvent polarization, and bath oscillators, respectively. Within our mixed quantum-classical approach, the proton and electron are treated quantum mechanically, whereas the solvent polarization and bath modes are treated classically. The protonic potential energy, Vp, is a double-well potential given by

dqpϕpi (qp)ϕpk (qp)∑ clα1 ′(0)cmα1(0) lm

⎤ e−βEα(Q ) 0 × dqpϕpl (qp)ϕpm(qp)⎥ −∞ ⎥⎦ (2π ℏ)ν ZQ

pp2

∫

(20)

where the equilibrium density of the PCET reactant state is defined as 1 nAeq = ∑ dX(|ϕe1⟩⟨ϕe1|h(−qp))αα ′ e−βHα (2π ℏ)ν ZQ αα ′

Vp(qp) = −

mpωp2 2

qp2 +

mp2ωp4 16V0

qp4 − λqp3

(23)

where ωp is the proton vibrational frequency and V0 is the energy difference between the minimum and maximum points of the potential. The potential energy matrix in the diabatic representation of the donor and acceptor electronic states is given by

∫

(21)

As can be seen, the resulting expression for kAB(t) has an equilibrium part (i.e., the terms that do not depend on time) and a dynamical part (i.e., the terms that depend on time). In the Simulation Details section, we explain how the sampling of the thermal equilibrium distributions is carried out. After sampling the initial conditions, the dynamical part is evaluated

⎛ V11(q ) + Vep(q ) ⎞ V12(qs) p p ⎜ ⎟ Ve(qp , qs) = ⎜ V22(qs) − Vep(qp)⎟⎠ ⎝ V12(qs) 3023

(24)


Article

Journal of Chemical Theory and Computation Table 1. Parameters for the Model Hamiltonian in Eq 22a

a

ms = 22 000

ωs = 3.72 × 10−4

qs1 = −2.13

qs2 = 2.13

mp = 1836.1 μ1 = 1.1 × 10−3 N = 12

ωp = 1.04 × 10−2 V12 = 2.45 × 10−2 Mj = ms

V0 = 1.2 × 10−2 ϕ = 8.0 η = msωs

μ2 = 5.84 × 10−3 λ = 0.0 T = 300

All values are in atomic units except temperature (T), which is in kelvin.

where the electron−proton coupling is given by Vep(qp) = μ2 tanh(ϕqp), the solvent polarization potentials in the donor and 1 acceptor states are given by V11(qs) = 2 msωs2(qs − qs )2 and

purposes of the calculation, the mean of qs is chosen such that the initial values of qs correspond to reactant configurations. Now that initial values of the positions and velocities have been selected, one needs to then sample from e−βEα(Q)/Z for each value of α. This can be done by first considering that we are interested in averages of the form:

1

1

V22(qs) = 2 msωs2(qs − qs )2 , respectively, and V12(qs) is the 2

diabatic coupling between the donor and acceptor states. The proton−solvent coupling is given by Vps(qp, qs) = −μ1qs tanh(ϕqp). The harmonic bath potential, Vb, is given by ⎛ cjqs ⎞ 1 ⎟ Vb = ∑ Mjωj2⎜⎜Q j − 2 j Mjωj2 ⎟⎠ ⎝

∑ ∫ dX ··· αα ′

where Z = ∑α∫ dX e . For convenience, we illustrate the sampling procedure for the case of a two-state system, although this can be generalized for any number of states. In this case, we may write the average explicitly as

(25)

∫ dX ··· e

(26)

+

where η denotes the strength of the coupling between the solvent polarization and bath modes, and ωc is the cutoff frequency. Discretizing this spectral density yields N harmonic oscillators with frequencies70 ⎛ j − 0.5 ⎞ ⎟ ωj = −ωc ln⎜ ⎝ N ⎠

−βE1(Q )

Z

∫ dX ··· e

+

∫ dX ··· e

−βE1(Q )

+

Z

∫ dX ··· e

−βE 2(Q )

Z

−βE 2(Q )

(30)

Z

where the first two terms correspond to αα′ = 11 and 12 and the latter two correspond to αα′ = 21 and 22. Let us first consider the 11 term. It may be rewritten as follows:

∫ dX ···

(27)

and coupling constants ⎛ 2ηMjωc ⎞1/2 cj = ωj⎜ ⎟ ⎝ Nπ ⎠

(29)

−βEα(Q)

where Mj and ωj are the mass and frequency of bath oscillator j, respectively. The bath is characterized by an ohmic spectral density: J(ω) = ηω exp( −ω/ωc)

e−βEα(Q ) Z

e−βE1(Q )

∫ dX(e−βE1(Q ) + e−βE2(Q ))

∫ dX ··· (28)

=

e−βE1(Q )

∫ dX e−βE1(Q )(1 + e−β(E2 − E1))

×

∫ dX e−βE1(Q ) ∫ dX e−βE1(Q )

⎛ ∫ dX ··· e−βE1(Q ) ⎞⎛ ∫ dX (1 + e−β(E2 − E1))e−βE1(Q ) ⎞−1 ⎟⎜ ⎟ = ⎜⎜ ⎟ −βE1(Q ) ⎟⎜ ∫ dX e−βE1(Q ) ⎝ ∫ dX e ⎠⎝ ⎠

In ref 66, the authors studied several parameter sets to mimic the different PCET mechanisms (i.e., concerted or sequential), which were based on the parameter sets in ref 38. In this study, we used the parameter set for Model I in ref 66 (see Table 1), which is meant to represent a symmetric hydrogen-bonding interface, for example, a dicarboxylic acid interface. Although this parameter set does not correspond to any particular system, it does yield some of the fundamental physics underlying this class of PCET reactions. In general, one may parametrize such model Hamiltonians for systems of interest by obtaining the parameters from electronic structure calculations and/or experiments.

(31)

Thus, the 11/12 term may be calculated as follows. One generates an ensemble of initial values of the positions and velocities. For each set of initial positions and velocities, one equilibrates the system on the ground-state surface, E1(Q), at the desired temperature using a Nosé−Hoover chain thermostat.71 A Nosé−Hoover chain thermostat is needed because the dynamics of the classical harmonic oscillators is nonergodic when a Nosé−Hoover thermostat is used. This generates an ensemble of initial conditions with the proper canonical distribution, from which microcanonical trajectories are run to propagate the time-dependent product species operator. The 11/12 term may be assembled by calculating the ground-state average within the first set of brackets in eq 31 and then reweighting it by the inverse of the ground-state average of (1 + e−β(E2−E1)) within the second set of brackets. One can follow a similar procedure for calculating the 21/22 term, where now the equilibration is performed on the E2(Q) surface and the first-excited-state average is reweighted by the inverse of the first-excited-state average of (1 + e−β(E1−E2)). As per eq 20, one must also calculate ∑lmcαl1′(0)cαm1(0)∫ 0−∞ dqpϕlp(qp)ϕmp (qp) at the beginning of each microcanonical trajectory.

4. SIMULATION DETAILS 4.1. Sampling the Thermal Equilibrium Distributions. To generate properly distributed initial conditions for calculating the time-dependent rate coefficient, one must sample the thermal equilibrium distributions, e−βEα(Q)/Z, for each value of α. (N.B.: Here, Q = {qs, Qj}.) We obtain initial values of the classical positions and velocities of the solvent polarization and bath oscillators from Gaussian distributions with standard deviations: σ v s = (βm s ) −1/2 and σ q s = −1/2 ω−1 for the solvent polarization and σVj = (βMj)−1/2 s (βms) −1/2 and σQj = ω−1 for the bath coordinates. For the j (βMj) 3024


Article

Journal of Chemical Theory and Computation It should be noted that the evaluation of the equilibrium density of the reactant species, nAeq, involves the same considerations as described in this subsection. 4.2. Simulating the Time Evolution of the Product Species Operator. After sampling the initial conditions, one must simulate the time evolution of the product species i k operator by evaluating ∑ikcαi2′(t)cαk2(t)∫ ∞ 0 dqp ϕp(qp)ϕp(qp) in eq 20. In particular, the time evolution of cαi2′(t) cαk2(t) is simulated using the sequential short-time propagation (SSTP) solution of the MQCL equation.60,61 This algorithm generates an ensemble of surface-hopping trajectories, containing segments in which the classical bath evolves either on single adiabatic surfaces or on the mean of two different surfaces [see Appendix 3 for the theoretical details of the SSTP evolution of cαi2′(t)cαk2(t)]. At each step of the SSTP algorithm, one needs to solve the timeindependent Schrödinger equation: h|̂ α ; qs⟩ = Eα(qs)|α ; qs⟩

(32)

where ĥ is the model Hamiltonian given in eq 22 without the purely classical terms. As mentioned earlier, this first involves expanding adiabatic wave function |α; qs⟩ in an orthonormal set of two-particle basis functions: |α ; qs⟩ =

α |ϕpm⟩|ϕen⟩ ∑ cmn

Figure 2. Ground-state and first-excited-state PESs (blue lines) and their corresponding mean PES (red line) as a function of qs for the PCET model (upper panel). The absolute value of the nonadiabatic coupling between these two surfaces, |d12|, as a function of qs (lower panel).

(33)

m,n

|ϕmp ⟩

where is chosen to be the mth eigenfunction of the quantum harmonic oscillator for a proton and |ϕne⟩ is the nth electronic diabatic state, where n varies from 1 (donor) to 2 (acceptor). Thus, for M protonic basis functions, the 2M × 2M Hamiltonian matrix is composed of matrix elements given by h{ij}{kl} = ⟨ϕpi ϕej|h|̂ ϕpk ϕel⟩, where the protonic integrals may be computed numerically using

using the SSTP algorithm. We first calculated the rate coefficient using adiabatic dynamics (i.e., d12 = 0) to investigate the effect of nonadiabatic transitions on the rate constant. This calculation involved 5 × 105 trajectories with a time step of 20 au (i.e., 0.48 fs). As can be seen in Figure 3, the rate coefficient

ϕpi (qp) = ⟨qp|ϕpi ⟩ 2

o

= (2ii! π )−1/2 αp1/2 e−αp (qp − qp )/2Hi(αp(qp − qpo)) (34) −1

Here, Hi(q) is the ith Hermite polynomial, αp = 3.97 Å , and qop = 0 Å. The integrals are calculated using the trapezoid rule with 300 points evenly spaced over the range −0.79 < qp < 0.79 Å. It was determined that M = 20 protonic states were sufficient for converging the eigenvalues. After constructing the Hamiltonian matrix, the eigenvalue problem, hc = Ec, is solved at each step of the simulation to obtain {cαmn} and {Eα}.

Figure 3. Rate coefficient as a function of time, kAB(t), calculated using adiabatic dynamics. The dotted red line indicates the plateau value.

5. RESULTS AND DISCUSSION In our calculation of kAB(t), we considered the participation of the two lowest adiabatic PESs. In the upper panel of Figure 2, we show the ground-state (11) and first-excited-state (22) PESs as a function of the solvent polarization coordinate, qs, for this model. We see that the reactant−product energy barrier on the ground PES is 1.0 kcal/mol (=1.68 kT). A minimum energy gap of 0.5 kcal/mol at qs = 0 Å separates the ground-state and first-excited-state PESs. The absolute value of the nonadiabatic coupling between these two surfaces, d12 = ⟨1; qs|∇qs|2; qs⟩, as a function of qs is depicted in the lower panel of Figure 2. As expected, it is peaked at qs = 0 Å, indicating that nonadiabatic transitions are most probable in the barrier-top region. Next, we present the results of the time-dependent rate coefficient, kAB(t), in eq 20 for the PCET model, calculated

first increases to its transition-state theory value, followed by a decrease to a plateau of kAB(t ∼ 9000 au) = 0 au (depicted by the red line). Thus, our adiabatic calculation yields a rate constant of 0 au, indicating that nonadiabatic transitions are essential for capturing the full magnitude of the rate constant. This behavior can be explained by considering the diagonal and off-diagonal contributions to the rate coefficient. When the system is initiated on a diagonal surface (i.e., α = α′) and remains on that surface, the imaginary part of the phase factor (which factors into the MQCL evolution of the product species operator) remains zero for all times. Therefore, there is no contribution to the rate coefficient from the trajectories that evolve on the 11 and 22 surfaces. On the other hand, when the system is initiated on an off-diagonal surface (i.e., α ≠ α′) and 3025


Article

Journal of Chemical Theory and Computation remains on that surface, the imaginary part of the phase factor will oscillate about zero. Therefore, the evolution on the 12 and 21 mean surfaces will give rise to a nonzero contribution to the rate coefficient initially (as seen in Figure 3). These contributions will have opposite signs because the imaginary part of the phase factor is a sine function and E1 − E2 < 0 and E2 − E1 > 0. However, because of the absence of an energy barrier on the 12 surface (see Figure 2), trajectories will readily oscillate between the reactant and product solvent configurations causing both the phase factor associated with each trajectory to evolve with different frequencies, ωαα′ = (Eα − Eα′)/ℏ, and the eigenvector coefficients in eq 20 to change rapidly. Together, this will cause the kAB(t) expectation value to decay to zero at longer times. When nonadiabatic transitions are taken into account, the situation changes. The corresponding time-dependent rate coefficient is shown in Figure 4. This result was calculated using

Figure 5. Contributions from the nonadiabatic trajectories starting on the 11, 12, 21, and 22 surfaces to the total time-dependent rate coefficient. The sum of the 12 and 21 contributions is also shown.

6. CONCLUSIONS We derived a novel formula for computing the rate constant of a PCET reaction, which starts from the time integral of a quantum mechanical flux−flux correlation function for a general reaction A ⇌ B.67−69 By introducing reactant and product species operators that are suitable for monitoring a PCET reaction, we obtained an expression for the timedependent rate coefficient. The resulting expression contains a dynamical part involving the time evolution of the product species operator and an equilibrium part. The dynamical part is approximated using MQCL dynamics, and the calculation of the full quantum equilibrium structure is simplified greatly by making a high-temperature approximation involving the classical DOF (which is a very good approximation because we are interested in simulating thermal PCET reactions). MQCL dynamics is used to evolve the product species operator because it treats decoherence more rigorously than do the FSSH-based methods. The final formula is general, requiring neither any prior assumptions based on the parameters of the system nor the mechanism of the reaction. Therefore, it may be applied to a wide range of vibronic couplings and solvent relaxation times, distinguishing it from golden-rule-based expressions, which rely on weak vibronic coupling and fast solvent relaxation. We then applied this PCET rate theory to a more complex version of the PCET model investigated in our previous studies,63,64 where now the solvent polarization coordinate is coupled to a dissipative bath of harmonic oscillators. Simulations of the time-dependent rate coefficient using both adiabatic and nonadiabatic dynamics revealed the importance of nonadiabatic dynamics starting on the mean PESs in accurately evaluating the rate constant. Our nonadiabtic rate constant was found to be in good agreement with the known exact result. These findings hold promise for accurately calculating PCET rate constants in more realistic reactions of chemical and biological interest and for determining which factors in the classical environment affect their rates.

Figure 4. Rate coefficient as a function of time, kAB(t), calculated using nonadiabatic dynamics. The dotted red line indicates the plateau value.

50 × 106 trajectories with a time step of 100 au (i.e., 2.4 fs) and an observable-cutting filter61 to reduce the statistical fluctuations. As can be seen, the nonadiabatic rate coefficient exhibits a similar behavior to the adiabatic one initially but eventually plateaus to a nonzero value of 1.34 × 10−6 au, yielding a rate constant of kAB = 1.34 × 10−6 au (i.e., 8.82 × 109 s−1). Comparing this value of the rate constant to the exact value of 1.71 × 10−6 au (i.e., 1.12 × 1010 s−1) reported in ref 66, we see that the results are in good agreement, attesting to the reliability and applicability of our PCET rate theory. There are several possible reasons for the discrepancy between our result and the exact result, including the approximations that are inherent in the MQCL approach and SSTP algorithm and the use of only two states in the calculation. To gain insight into the contributions from the trajectories starting on the various surfaces to the time-dependent nonadiabatic rate coefficient, we decomposed it into its 11, 12, 21, and 22 components. As can be seen in Figure 5, the rate coefficient takes on a nonzero value, essentially because of nonadiabatic dynamics starting on the 12 and 21 mean surfaces. The 11 component is comparably very small and is responsible for the noise observed in the total rate coefficient, whereas the 22 component is zero. Thus, in this case, it can be said that the total rate constant is the sum of the 12 and 21 components (also shown in Figure 5). These results are consistent with our findings in refs 63 and 64, which demonstrated the important role of the mean PESs in the PCET reaction dynamics of the models studied.

■

APPENDIX 1. EVALUATION OF MATRIX ELEMENTS IN WαA′α(X, 0) ̂

Matrix element ⟨2Q ′ − Q − (Z /2)|e−βH |Q − Z /2⟩ appearing in eq 12 may be evaluated as follows. The Hamiltonian of the total system, Ĥ , may be partitioned into the kinetic energy of the environment and the remainder, ĥ. On the basis of this ̂ partitioning, one may factorize e−βH using a symmetric Trotter 72 decomposition as 3026


Article

Journal of Chemical Theory and Computation ̂

̂

̂

2

e−βH = e−β(h /2) e−β ∑i(Pi /2Mi) e−β(h /2) + 6(β 3)

β /2Mi [Pi + (iMi /β ℏ)(2Q i − 2Q i′)]

ui =

(35)

Inserting this into the matrix element yields Z −βĤ Z |e |Q − ⟩ 2 2 Z −β(h ̂/2) −β ∑i(Pi2 /2Mi) −β(h ̂/2) Z = ⟨2Q ′ − Q − |e |Q − ⟩ e e 2 2 Z −β(h ̂/2) −β ∑i (Pi2 /2Mi) = dQ 1 dQ 2⟨2Q ′ − Q − |e |Q 1⟩⟨Q 1|e |Q 2⟩ 2 ̂ Z × ⟨Q 2|e−β(h /2)|Q − ⟩ 2 Z −β(h ̂/2)(2Q ′− Q − Z /2) = dQ 1 dQ 2 e ⟨2Q ′ − Q − |Q 1⟩ 2

dui =

β /2Mi dPi , the integrals over Pi may be solved to

obtain

2Miπ /β . Hence, the above expression simplifies to

⟨2Q ′ − Q −

=

APPENDIX 2. PERFORMING THE Z INTEGRATION IN WαA′α(X, 0) Using the fact that ∫ dZδ(Z)f(Z) = f(0), we may carry out the integration over Z in eq 18 to obtain

̂

2

̂

× e−β(h /2)(Q − Z /2)

Z ⟩ 2 (36)

W Aα ′ α(X , 0) =

1 (2π ℏ)ν ZQ

2

We next consider the ⟨2Q′ − Q − Z/2|e−β∑i(Pi /2Mi)|Q − Z/2⟩ term in eq 14. Rewriting it explicitly in terms of the product of states over all the classical DOFs, we obtain ⟨2Q ′ − Q −

Z |e 2

|Q −

∏ ⟨2Q i′ − Q i −

=

i

i

∏

= e

i

1 2π ℏ

⟨α1′; Q |α2 ; Q ⟩⟨α2 ; Q |α3; Q ⟩⟨α3; Q |α ; Q ⟩

W Aα ′ α(X , 0) =

∑ ∑

2 i

i

Zi |Pi⟩ 2

∑ ∑ =

i

∏

= e =

1 2π ℏ

2 i

i

i

i

=

∫ dPi

=

■

2

i

i

i

(37)

To arrive at the final line of eq 37, we inserted a complete set of momentum states, ∫ dPi|Pi⟩⟨Pi|, and then used the fact that ⟨Q |P⟩ = (1/ 2π ℏ )e(i / ℏ)P·Q . By letting

+∞

−∞

∑ ciα1 ′ckα1∫ i,k

∫ dPie−(β/2M )(P+(iM /βℏ)(2Q −2Q ′)) i

∑ ciα1 ′ckα1∫ i,k

1 −(Mi /2βℏ2)(2Q i − 2Q i′)2 e 2π ℏ i

(41)

i,k

i −[(β /2Mi)(Pi + (iMi / β ℏ)(2Q i − 2Q i′))2 + (β /2Mi)(Mi2 / β 2ℏ2)(2Q i − 2Q i′)2 ]

∏

cijα ′cklα⟨ϕej|ϕe1⟩⟨ϕpi |h( −qp)|ϕpk ⟩⟨ϕe1|ϕel⟩

∑ ciα1 ′ckα1⟨ϕpi|h(−qp)|ϕpk⟩

∫ dPi e−(β/2M )[P +(2iM P/βℏ)(2Q −2Q ′)] i i

∑ ∑

Because of the two Dirac delta functions, this expression has a nonzero value if and only if j = l = 1. Therefore, we obtain

∫

i

cijα ′cklα⟨ϕej|⟨ϕpi ||ϕe1⟩⟨ϕe1|h( −qp)|ϕpk ⟩|ϕel⟩

i , k j , l = 1,2

e−(i/ ℏ)Pi(Q i − Zi /2) 2 1 =∏ dPi e−β(Pi /2Mi) e−(i/ ℏ)Pi(2Q i − 2Q i′) 2π ℏ i 1 2π ℏ

(40)

i , k j , l = 1,2

∫

∏

cijα ′cklα⟨ϕej|⟨ϕpi ||ϕe1⟩⟨ϕe1|h( −qp)|ϕpk ⟩|ϕel⟩

We may now simplify the sums by rewriting them as

e−(i/ ℏ)Pi(Q i − Zi /2) 2 1 =∏ dPi e−β(Pi /2Mi) e(i/ ℏ)Pi(2Q i′− Q i − Zi /2) i 2π ℏ

=

e−βEα(Q ) (2π ℏ)ν ZQ

i , k j , l = 1,2

∫ dPi e−β(P /2M )⟨2Q i′ − Q i −

(39)

The four Dirac delta functions in this equation tell us that nonzero contributions arise when α = α3 = α2 = α1′ and α1 = α′. Finally, WαA′α(X, 0) simplifies to

Zi −β(Pi2 /2Mi) |e |Pi⟩ 2

∫ dPi⟨2Q i′ − Q i −

cijα1cklα1⟨ϕej|⟨ϕpi ||ϕe1⟩⟨ϕe1|h( −qp)|ϕpk ⟩|ϕel⟩

× e−(β /2)Eα2(Q ) e−(β /2)Eα3(Q ) × ⟨α′; Q |α1; Q ⟩

Zi −β(Pi2 /2Mi) Z |e |Pi⟩⟨Pi|Q i − i ⟩ 2 2

i −(i/ ℏ)Pi(Q i − Zi /2)

∏

=

1 2π ℏ

α1α1′α2α3 i , k j , l = 1,2

Zi −β(Pi2 /2Mi) Z |e |Q i − i ⟩ 2 2

∏ ∫ dPi⟨2Q i′ − Q i −

=

′

∑ ∑ ∑

Z ⟩ 2

−β ∑i (Pi2 /2Mi)

(38)

■

× ⟨Q 1|e−β ∑i(Pi /2Mi)|Q 2⟩⟨Q 2|Q − Z /2⟩e−β(h /2)(Q − Z /2) ̂

∏ i

∫

= e−β(h /2)(2Q ′− Q − Z /2)⟨2Q ′ − Q − Z /2|e−β ∑i(Pi /2Mi)|Q −

Z −β ∑i(Pi2 /2Mi) Z |e |Q − ⟩ 2 2 2 2 Mi e−(Mi /2βℏ )(2Q i − 2Q i′) 2 2βπ ℏ

⟨2Q ′ − Q −

∫

2

and

0

−∞

dqpϕpi (qp)h( −qp)ϕpk (qp)

dqpϕpi (qp)ϕpk (qp)

(42)

APPENDIX 3. THEORETICAL DETAILS OF THE MQCL TIME EVOLUTION OF CαI2′(T)CαK2(T)

The time evolution of cαi2′(t)cαk2(t) ≡ Cαα ik ′(X, t) is carried out using the SSTP solution of the MQCL equation, as dictated by60,61 3027


Article

Journal of Chemical Theory and Computation L

Cikαα ′(X ,

∑

t) ≈

(α1α1′) ...(αLαL′)

(9) Warren, J.; Tronic, T.; Mayer, J. Chem. Rev. 2010, 110, 6961. (10) Costentin, C. Chem. Rev. 2008, 108, 2145. (11) Song, N.; Gagliardi, C. J.; Binstead, R. A.; Zhang, M.-T.; Thorp, H.; Meyer, T. J. J. Am. Chem. Soc. 2012, 134, 18538. (12) Roubelakis, M. M.; Bediako, D. K.; Dogutan, D. K.; Nocera, D. G. Energy Environ. Sci. 2012, 5, 7737. (13) Schrauben, J. N.; Hayoun, R.; Valdez, C. N.; Braten, M.; Fridley, L.; Mayer, J. M. Science 2012, 336, 1298. (14) Pizano, A. A.; Yang, J. L.; Nocera, D. G. Chem. Sci. 2012, 3, 2457. (15) Cukier, R. J. Phys. Chem. 1996, 100, 15428. (16) Cukier, R.; Nocera, D. Annu. Rev. Phys. Chem. 1998, 49, 337. (17) Cukier, R. J. Phys. Chem. A 1999, 103, 5989. (18) Soudackov, A.; Hammes-Schiffer, S. J. Chem. Phys. 2000, 113, 2385. (19) Hammes-Schiffer, S.; Soudackov, A. J. Chem. Phys. 2008, 112, 14108. (20) Siegbahn, P. E.; Blomberg, M. Chem. Rev. 2010, 110, 7040. (21) Kretchmer, J.; Miller, T. F., III J. Chem. Phys. 2013, 138, 134109. (22) Migliore, A.; Polizzi, N. F.; Therien, M. J.; Beratan, D. N. Chem. Rev. 2014, 114, 3381−3465. (23) Marcus, R. J. Chem. Phys 1956, 24, 966. (24) Marcus, R. Can. J. Chem. 1959, 37, 155. (25) Borgis, D.; Lee, S.; Hynes, J. T. Chem. Phys. Lett. 1989, 162, 19. (26) Borgis, D.; Hynes, J. T. J. Chem. Phys. 1991, 94, 3619. (27) Kuznetsov, A. M.; Ulstrup, J. Can. J. Chem. 1999, 77, 1085. (28) Cukier, R. J. Phys. Chem. 1994, 98, 2377. (29) Georgievskii, Y.; Stuchebrukhov, A. A. J. Chem. Phys. 2000, 113, 10438. (30) Soudackov, A.; Hammes-Schiffer, S. J. Chem. Phys. 1999, 111, 4672. (31) Soudackov, A.; Hatcher, E.; Hammes-Schiffer, S. J. Chem. Phys. 2005, 122, 014505. (32) Hatcher, E.; Soudackov, A.; Hammes-Schiffer, S. Chem. Phys. 2005, 319, 93. (33) Hammes-Schiffer, S.; Soudackov, A. J. Phys. Chem. B 2008, 112, 14108. (34) Soudackov, A.; Hammes-Schiffer, S. J. Chem. Phys. 2015, 143, 194101. (35) Tully, J.; Preston, R. J. Chem. Phys. 1990, 93, 1061. (36) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2004, 121, 3368. (37) Habershon, S.; Manolopoulos, D. E.; Markland, T. E.; Miller, T. F. Annu. Rev. Phys. Chem. 2013, 64, 387. (38) Fang, J.-Y.; Hammes-Schiffer, S. J. Chem. Phys. 1997, 106, 8442. (39) Fang, J.-Y.; Hammes-Schiffer, S. J. Chem. Phys. 1997, 107, 8933. (40) Fang, J.-Y.; Hammes-Schiffer, S. J. Chem. Phys. 1997, 107, 5727. (41) Kretchmer, J.; Miller, T. F., III Inorg. Chem. 2016, 55, 1022. (42) Miller, T. F., III J. Chem. Phys. 2008, 129, 194502. (43) Menzeleev, A. R.; Ananth, N.; Miller, T. F., III J. Chem. Phys. 2011, 135, 074106. (44) Menzeleev, A. R.; Bell, F.; Miller, T. F., III J. Chem. Phys. 2014, 140, 064103. (45) Subotnik, J. E.; Ouyang, W.; Landry, B. R. J. Chem. Phys. 2013, 139, 214107. (46) Landry, B. R.; Subotnik, J. E. J. Chem. Phys. 2011, 135, 191101. (47) Falk, M.; Landry, B.; Subotnik, J. J. Phys. Chem. B 2014, 118, 8108. (48) Kelly, A.; Markland, T. E. J. Chem. Phys 2013, 139, 014104. (49) Schwartz, B. J.; Bittner, E. R.; Prezhdo, O. V.; Rossky, P. J. J. Chem. Phys. 1996, 104, 5942. (50) Jasper, A. W.; Stechmann, S. N.; Truhlar, D. G. J. Chem. Phys. 2002, 116, 5424. (51) Bedard-Hearn, M. J.; Larsen, R. E.; Schwartz, B. J. J. Chem. Phys. 2005, 123, 234106. (52) Landry, B. R.; Subotnik, J. E. J. Chem. Phys. 2012, 137, 22A513. (53) Martens, C. C.; Fang, J. J. Chem. Phys. 1997, 106, 4918. (54) Donoso, A.; Martens, C. C. J. Chem. Phys. 1998, 102, 4291. (55) Donoso, A.; Martens, C. C. J. Chem. Phys. 2000, 112, 3980.

[∏ >αj−1αj′−1(t j − 1 , t j)e iLαj−1αj′−1Δt j j=1 ′

× (δαj−1αjδ αj′−1αj′ − Δt 1αj−1αj′−1, αjαj′)]CikαLαL(X )

(43)

where time interval t is divided into L segments such that the jth segment has length Δt j = t j − t j − 1 = Δt, >αj−1α′j−1(t j − 1 , t j) = e(i/ ℏ)(Eαj−1 − Eαj′−1)Δt j is the phase factor associated with the jth segment, the action of classical Liouville operator L is given by iLαα ′ =

P ∂ 1 ∂ · + (F α + F α ′) · M ∂Q 2 ∂P

(44)

and the action of nonadiabatic transition operator 1 is given by ⎛ ∂ ⎞ P 1 ·dαβ ⎜1 + Sαβ · ⎟δα ′ β ′ ⎝ ∂P ⎠ M 2 ⎛ ∂ ⎞ P 1 − ·dα*′ β ′⎜1 + Sα*′ β ′· ⎟δαβ ⎝ ∂P ⎠ M 2

1αα ′ , ββ ′ = −

(45)

In eqs 44 and 45, F α = −⟨α ; Q |∂h(̂ Q )/∂Q |α ; Q ⟩ is the Hellmann−Feynman force corresponding to PES α, dαβ = ⟨α ; Q |∂/∂Q |β ; Q ⟩ is the nonadiabatic coupling matrix element, and Sαβ = (Eα − Eβ)dαβ(P/M·dαβ)−1. Thus, the classical DOFs are propagated using the Hellmann−Feynman forces either on the PES corresponding to one of the states, α, with energy Eα or on the mean PES corresponding to states α and α′ with energy (Eα + Eα′)/2. In practice, the momentumjump approximation61,73 is used to compute the action of 1αα ′ , ββ ′, which leads to changes in the momenta of the classical environment caused by the energy exchange between the quantum and classical subsystems during nonadiabatic transitions. The algorithmic details of implementing this solution are given in refs 60 and 61.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by grants from the Natural Sciences and Engineering Research Council of Canada and Alberta Innovates Technology Futures.

■

REFERENCES

(1) Huynh, M.; Meyer, T. Chem. Rev. 2007, 107, 5004. (2) Meyer, T.; Huynh, M.; Thorp, H. H. Angew. Chem., Int. Ed. 2007, 46, 5284. (3) Hammes-Schiffer, S.; Stuchebrukhov, A. Chem. Rev. 2010, 110, 6939. (4) Weinberg, D.; Gagliardi, C.; Hull, J.; Murphy, C. F.; Kent, C.; Westlake, B.; Paul, A.; Ess, D.; McCafferty, D. G.; Meyer, T. Chem. Rev. 2012, 112, 4016. (5) Kaila, V. R. I.; Verkhovsky, M. I.; Wikstrom, M. Chem. Rev. 2010, 110, 7062. (6) Babcock, G. T.; Wikström, M. Nature 1992, 356, 301. (7) Milligan, J. R.; Aguilera, J. A.; Hoang, O.; Ly, A.; Tran, N. Q.; Ward, J. F. J. Am. Chem. Soc. 2004, 126, 1682−1687. (8) Petek, H.; Zhao, J. Chem. Rev. 2010, 110, 7082. 3028


Article

Journal of Chemical Theory and Computation (56) Santer, M.; Manthe, U.; Stock, G. J. Chem. Phys. 2001, 114, 2001. (57) Kapral, R.; Ciccotti, G. J. Chem. Phys. 1999, 110, 8919. (58) Nielsen, S.; Kapral, R.; Ciccotti, G. J. Stat. Phys. 2000, 101, 225. (59) Wan, C.; Schofield, J. J. Chem. Phys. 2000, 113, 7047. (60) MacKernan, D.; Kapral, R.; Ciccotti, G. J. Phys.: Condens. Matter 2002, 14, 9069. (61) Hanna, G.; Kapral, R. J. Chem. Phys. 2005, 122, 244505. (62) MacKernan, D.; Ciccotti, G.; Kapral, R. J. Phys. Chem. B 2008, 112, 424. (63) Shakib, F.; Hanna, G. J. Chem. Phys. 2014, 141, 044122. (64) Shakib, F.; Hanna, G. J. Chem. Phys. 2016, 144, 024110. (65) Yamamoto, T. J. Chem. Phys. 1960, 33, 281. (66) Ananth, N.; Miller, T. F., III Mol. Phys. 2012, 110, 1009. (67) Kim, H.; Kapral, R. J. Chem. Phys. 2005, 123, 194108. (68) Hanna, G.; Kim, H.; Kapral, R. In Quantum Dynamics of Complex Molecular Systems; Micha, D. A., Burghardt, I., Eds.; Springer: Berlin, 2006. (69) Kim, H.; Kapral, R. J. Chem. Phys. 2005, 122, 214105. (70) Craig, I. R.; Manolopoulos, D. E. J. Chem. Phys. 2005, 122, 084106. (71) Tuckerman, M. Statistical Mechanics: Theory and Molecular Simulation; Oxford University Press: Oxford, 2010. (72) Trotter, H. F. Proc. Am. Math. Soc. 1959, 10, 545. (73) Sergi, A.; MacKernan, D.; Ciccotti, G.; Kapral, R. Theor. Chem. Acc. 2003, 110, 49−58.

3029


Mixed Quantum-Classical Liouville Approach for Calculating Proton

Recommend Documents