Mixed Quantum-Classical Study of Nonadiabatic Curve Crossing in

Subscriber access provided by The University of British Columbia Library

Article

Mixed Quantum-Classical Study of NonAdiabatic Curve Crossing in Condensed Phases Weiwei Xie, Meng Xu, Shuming Bai, and Qiang Shi J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b11695 • Publication Date (Web): 03 Feb 2016 Downloaded from http://pubs.acs.org on February 11, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Mixed Quantum-Classical Study of Non-Adiabatic Curve Crossing in Condensed Phases Weiwei Xie, Meng Xu, Shuming Bai, and Qiang Shi∗ Beijing National Laboratory for Molecular Sciences, State Key Laboratory for Structural Chemistry of Unstable and Stable Species, Institute of Chemistry, Chinese Academy of Sciences, Zhongguancun, Beijing 100190, China E-mail: [email protected] Phone: +86 10 82616163. Fax: +86 10 82616163

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We apply the mixed quantum-classical Liouville (MQCL) equation to investigate the non-adiabatic curve crossing in condensed phases. More specifically, electron transfer rate constants of the spin-Boson model are calculated by employing a rate constant expression using the collective solvent polarization as the reaction coordinate. In the calculation, classical nuclear degrees of freedom are initially sampled at the transition state configuration, and initial state for the electronic degree of freedom is obtained from a mixed quantum-classical Boltzmann distribution. Different contributions to the electron transfer rate from the diagonal and off-diagonal elements of the initial density matrix, and contributions from trajectories with positive and negative initial velocities are analyzed. It is shown that the off-diagonal elements of the initial density matrix play an important role in the total electron transfer rate. The MQCL results are also compared with those calculated using Ehrenfest dynamics. It is found that, although the Ehrenfest dynamics is inaccurate when using the reactive flux rate expression directly, it can give reasonably accurate results when calculating individual contributions from the diagonal and off-diagonal elements of the initial density matrix.

Keywords mixed quantum-classical Liouville equation, reactive flux method, Ehrenfest dynamics, spinBoson

1

Introduction

The phenomena of non-adiabatic curve crossing are ubiquitous in condensed phase chemical processes, such as charge and energy transfer reactions, 1 and photochemical reactions. 2 Unlike the case of adiabatic reactions where classical dynamics can be applied if the quantum effect of the nuclear degrees of freedom can be neglected, classical dynamics is not adequate in describing the coupled electron-nuclear motion in non-adiabatic reactions, while satisfactory 2


Page 2 of 30

Page 3 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mixed quantum-classical descriptions of non-adiabatic dynamics still present a significant challenge in theoretical chemistry. Many existing theories can be applied to non-adiabatic reactions, such as the Marcus theory 3,4 and the more general Fermi’s golden rule (FGR). 5 But they are based on pertubative treatment of the electronic coupling between different electronic states. To go beyond the purturbative treatment of the electronic coupling, the Landau-Zener formula of curvecrossing 6,7 has been widely used. For example, Cline and Wolynes 8 proposed a semiclassical description of non-adiabatic curve-crossing using surface hopping and the Landau-Zener formula. Straub and Berne 9 have derived rate expressions that further include the effect of solvent relaxation. More recently, Zhao et al. 10,11 have derived non-adiabatic rate expressions using the more accurate Zhu-Nakamura theory 12 instead of the Landau-Zener formula. Mixed quantum-classical simulations have also been used to calculate non-adiabatic reaction rate constants without invoking either the Landau-Zener formula or the Zhu-Nakamura theory. Such simulations can be classified into two categories. The first class of mixed quantum-classical methods are based on direct simulation of the non-adiabatic dynamics with the nuclear configurations sampled from the equilibrated donor state. By using such direct simulation method, Subotnik and coworkers have analyzed in detail how the popular fewest switch surface hopping (FSSH) method should be modified to obtain the accurate electron transfer (ET) reaction rate constants. 13–15 Xie et al. have found that many approximated methods including the Ehrenfest dynamics and surface hopping methods can give reasonable ET reaction rate constants in a wide range of parameter regimes. 16 Another class of methods to calculate non-adiabatic rate constants are based on the reactive flux method, 17–19 in which the initial nuclear coordinates are sampled at the transition state configuration. In the case of adiabatic reactions, the reactive flux method successfully overcomes the rare event sampling problem, and connects naturally to the celebrated transition state theory. 18 The reactive flux method in combination with the mixed quantumclassical simulations is also widely applied in calculating non-adiabatic reaction rate con-

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

stants. For example, Hammes-Schiffer and Tully have developed a method to calculate nonadiabatic reaction rate constants using the molecular dynamics with quantum transitions (MDQT) based on surface hopping. 20 This method was later used by Kim and HammesSchiffer to calculate the proton transfer (PT) rate of a double well model coupled to harmonic bath, 21 which is compared with benchmark calculations using the quasi-adiabatic propagator path integral (QUAPI) method. 22 More recently, Jain and Subotnik 23 calculated rate constants of the spin-Boson model using a similar method based on the augmented fewest switch surface hopping (A-FSSH) method, and have shown the importance of properly treating decoherence and forbidden hops. Kapral and coworkers have combined the non-adiabatic reactive flux method and the mixed quantum-classical Liouville (MQCL) equation in simulating reaction dynamics in the ET 24–26 and PT 27–30 reactions. Ka and Thompson 31 have also studied PT in the Azzouz-Borgis (AB) model 32 by employing the Meyer-Miller mapping method 33 and the non-adiabatic reactive flux formalism. In this work, we calculate the ET rate constants of the spin-Boson model using the MQCL method and the reactive flux method. In difference to previous MQCL studies of non-adiabatic rate dynamics by Kapral and coworkers, 24–26 the MQCL equation in the diabatic basis is employed in this work, and the collective polarization coordinate is used as the reaction coordinate. As the MQCL equation is exact for the spin-Boson model, an important goal of the current study is to reveal contributions from different terms of the initial density matrix at the transition state, and to use the simulation result as benchmarks comparing with results from the approximate method using Ehrenfest dynamics. The remaining sections of this paper are arranged as follows: in section 2, we present the model Hamiltonian, the reactive flux method to calculate non-adiabatic reaction rate constants, and the MQCL equation in the diabatic basis. In section 3, we perform MQCL simulations to obtain the non-adiabatic ET rate constants, and compare the results with those from the Ehrenfest dynamics. Conclusions and discussions are made in section 4.

4


Page 4 of 30

Page 5 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2 2.1

Theory Model Hamiltonian

As in several recent studies, 13–16,34 we employ the spin-Boson Hamiltonian 35,36 as a microscopic model of ET reactions. The total Hamiltonian describes a two level system coupled linearly to a bath of harmonic oscillators, which is written as " 2 # N X p2j ǫ 1 c σ j z H = V σx + σz + , + mj ωj2 xj + 2 2m 2 mj ωj2 j j=1

(1)

where σx and σz are the Pauli matrices, σx = |diha| + |aihd| and σz = |dihd| − |aiha|; |di and |ai denote the donor and acceptor states, V is the electronic coupling, and ǫ is the energy difference between them; xj , pj , and ωj are the coordinate, momentum, and frequency of the jth bath mode; the harmonic bath couples linearly to the electronic degree of freedom and causes energy fluctuations of the |di and |ai states, where cj is the coupling coefficient between the jth bath coordinate and the system operator σz . P The collective coordinate X is defined as X = ci xi , which will be used as the reaction

coordinate of the ET reaction in the derivations below. The system-bath coupling in the spin-Boson model system is characterized by the spectral density J(ω) defined as

J (ω) =

2.2

π X c2j δ(ω − ωj ) . 2 j mj ωj

(2)

The reactive flux method

In this subsection, we briefly present the reactive flux method in calculating the ET rate constants of the spin-Boson model. We first assume that the ET reaction occurs in the normal regime, so the collective coordinate X = 21 ∆E = 21 (Hd − Ha − ǫ) is a good reaction coordinate, while possible extensions to inverted region will be discussed later. Here, Hd = P p2j P p2j 1 ǫ 1 ǫ 2 2 2 2 + + m ω x + c x , and H = − + + m ω x − c x j j a j j are the nuclear j 2mj j 2mj 2 2 j j j 2 2 j j j 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 30

Hamiltonian on the donor and acceptor states, respectively. In realistic systems, the above defined X corresponds to the collective solvent polarization. 19 By using the collective coordinate X as the reaction coordinate, the population operator for the product state can be defined as

h = Θ(X − X ∗ ) = Θ(X − X ∗ )(|dihd| + |aiha|) ,

(3)

where Θ is the Heaviside function, and X ∗ defines the transition state where HD = HA . It is noted that the |dihd| + |aiha| term in the above equation simply comes from the identity operator for the electronic degree of freedom. The product projection operator h defined in Eq. (3) is the same as those in Refs., 20,23 and has the advantage that the flux operator (see below) has a δ-function part that confines the sampling of nuclear degrees of freedom at the transition state. Although an alternative definition of the product projection operator has also been used in literature, 24 the corresponding flux operator will have a term that is not confined to the transition state configuration. We thus limit our discussions below to the projection operator defined in Eq. (3). We start from the linear response theory in deriving the reactive flux formalism. 37 A small perturbation proportional to the h operator is first added to the total equilibrium state, such that ρ(0) = e−β(H+f h) /Tre−β(H+f h) ,

(4)

where f is a small parameter. By treating the nuclear degrees of freedom classically and performing first order expansion of the small parameter f , the initial density matrix can be approximated as e−βH (1 − f βδh) , ρ(0) ≈ Z

(5)

where Z is the partition function of the unperturbed Hamiltonian. The deviation of the

6


Page 7 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


product population from its equilibrium value at time t is then given by

δP (t) ∝ −f βhδh(0)δh(t)i ,

(6)

where the average is taken over the Boltzmann distribution of the mixed quantum-classical system. The phenomenological (total) rate constant k is defined by using the relation d δP (t) dt

= −kδP (t), and can be calculated as 17–19

k = lim k(t) , t→τp

k(t) =

hF (0)δh(t)i , hδh(0)δh(0)i

(7)

where τp is a short interval for k(t) to reach a plateau, and the flux operator F is defined as X i F ≡ h˙ = [H, h] = δ(X − X ∗ ) cj vj , ~ j

(8)

where vj is the velocity of the jth bath mode. The time correlation function hF δh(t)i in Eq. (7) can be further written as Z∗ hF δh(t)i = Z

Z

dp

Z

dxρ0 (p, x)

X

cj vj δh(t) ,

(9)

j

h i P where Z ∗ = Tr e−βH δ( j cj xj − X ∗ ) is the partition function at the transition state, and P ∗ e−βHB δ( j cj xj − X ∗ ) e−β(HS +X σz ) h i , ⊗ ρ0 = TrS [e−β(HS +X ∗ σz ) ] Tr e−βHB δ(P c x − X ∗ ) B j j j

(10)

with HS = V σx + 2ǫ σz . The δ-function in the flux operator F plays an important role in calculating the rate constant for adiabatic (single surface) reactions, 17,18 as it allows to start the simulation by placing the nuclear coordinates at the transition state configuration, and thus effectively avoids the problem of rare event sampling. However, calculation of the transmission coeffi-

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

cients is slightly different for non-adiabatic reactions when coupled electron-nuclear motion on multiple surfaces is involved. Even with the initial nuclear configurations sampled at the transition state, the transmission coefficient may still be small (∝ V 2 ) in the non-adiabatic limit. To better understand this behavior in the non-adiabatic limit, we notice that the initial density matrix in Eq. (10) contains terms like |dihd| and |aiha| (diagonal terms) and |diha|, |aihd| (off-diagonal terms). For the |dihd| term which indicates that the initial electronic population is on the donor state, most of the trajectories will not reach the product state when V is small, even if the initial collective velocity X˙ > 0, so the transmission coefficient from this term is small. For the |aiha| term, most of the trajectories will reach the product state, no matter whether the initial collective velocity is positive or negative, so transmission coefficient from this term is also small due to the cancellation of contributions from trajectories with positive and negative initial collective velocities. Contributions from the |diha| and |aihd| terms are scarcely discussed in the literature. In this work, we will show that contributions from these terms are actually comparable with those from the diagonal terms. So even though the reactive flux formalism helps to resolve the problem of rare event sampling by placing the nuclear configuration at the transition state, another important problem is to accurately calculate the transmission coefficient through the non-adiabatic dynamics simulation.

2.3

MQCL equation in the diabatic basis

The MQCL equation has recently developed into an important method in simulating nonadiabatic dynamics in complex molecular systems. 38–47 As stated above, Kaparal and coworkers have applied the MQCL method and the reactive flux formalism in the adiabatic basis in simulating non-adiabatic reactions. 24–30 In this work, we employ the MQCL equations in the

8


Page 8 of 30

Page 9 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


diabatic basis, which can be obtained from the MQCL equation in the operator form: 39,41,42 ∂ρw i 1 = − [Hw , ρw ] + [{Hw , ρw } − {ρw , Hw }] , ∂t ~ 2

(11)

where ρw is the partially Wigner transformed density matrix, Hw is the partial Wigner transform of the total Hamiltonian H. For simplicity, we will drop the subscript w in the following derivations. For the spin-Boson model, the MQCL equation in the diabatic basis is given by 48 X ∂ρdd ∂ρdd ∂ρdd 2 , = −i∆(ρad − ρda ) − pj − (ωj xj + cj ) ∂t ∂x ∂p j j j X ∂ρaa ∂ρaa ∂ρaa 2 , = i∆(ρad − ρda ) − pj − (ωj xj − cj ) ∂t ∂xj ∂pj j ! X X ∂ρda ∂ρda ∂ρda 2 − ωj xj = −i ǫ + 2 cj xj ρda − i∆(ρaa − ρdd ) − pj ∂t ∂x ∂pj j j j ! X X ∂ρad ∂ρad ∂ρad 2 . = i ǫ+2 cj xj ρad − i∆(ρdd − ρaa ) − pj − ωj xj ∂t ∂x ∂p j j j j

(12a) (12b) , (12c) (12d)

To calculate the correlation function hF δh(t)i in Eq. (9), we use ρ0 defined in Eq. (10) as the initial density matrix, and calculate the correlation between the collective velocity X˙ at time zero and product population δh(t). The initial state in Eq.(10) have both diagonal and off-diagonal terms: ρ0 = ρ0dd |dihd| + ρ0da |diha| + ρ0ad |aihd| + ρ0aa |aiha|. We note that, the MQCL equation consists of a set of linear equations in the sense that, the initial density matrix can be decomposed into different terms, and their individual contributions can be added up to obtain the final result. In the MQCL calculations, the spectral density in Eq. (2) is decomposed into discreet harmonic oscillator modes. The initial positions and velocities of each harmonic oscillator mode are sampled at the transition state. The Z ∗ /Z term in Eq. (9) can be calculated using the method described in Appendix A, and the sampling of initial positions and velocities are performed using the method presented in Appendix B. We also use the trajectory based 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 30

Monte Carlo method in previous works 42,45 to propagate the MQCL equation. Such numerical methods are often less effective in longer simulation times, 44–46,49 but we found that it gives reasonably converged results in simulating the systems presented in the next section 3. Correlation function between the initial collective velocity and the product population h(t) is then calculated. In some of the calculations, we also use the electronic population operator on the acceptor state pa ≡ |aiha| instead of h(t) in Eq. (3) to define the product population. We find that this alternative definition of the product population operator works well in calculating the ET rate constants. Another advantage of this substitution is that it can be generalized to ET reactions in the inverted region.

3

Results

We first present numerical simulations using the MQCL method to investigate contributions from different terms of the initial density matrix to the total non-adiabatic ET rate. The Brownian oscillator spectral density is used, where ηω 1 . J(ω) = λΩ2 2 2 2 M (ω − Ω2 )2 + η 2 ω 2

(13)

Here M and Ω denote the mass and the frequency of the Brownian oscillator, λ is the reorganization energy, and η is the friction coefficient. The parameters are taken from two recent studies, 13,14 with the reorganization energy λ = 2.39×10−2 a.u. (atomic unit), T =300 K, except that the harmonic oscillator frequency Ω = 1.0 × 10−3 a.u which is slightly larger than that used in Refs., 13,14 the coupling strength V = 2.0 × 10−4 a.u., the energy difference ǫ = 0, and the friction constant η = 1.2 × 10−3 a.u. The time step is set to 12 a.u, and averages over 107 trajectories are used in the MQCL simulation. The number of harmonic bath modes to decompose the Brownian oscillator spectral density is 1000. Figure 1 shows the ground and excited state free energy profiles as a function of the reaction coordinate X. Equations to calculate the free energy profiles are given in Appendix 10


Page 11 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A. Figure 2 shows k(t) in Eq. (7) calculated by the MQCL method in the diabatic basis, from which the ET rate constant can be obtained. Contributions from different terms in the initial density matrix are also shown in Figure 2. The (total) ET rate constant calculated by FGR 1,36 using the same parameters is 6.34 × 10−8 a.u., which agrees well with the MQCL result. It can also be seen from Figure 2 that, contributions from the off-diagonal |diha| and |aihd| terms are comparable to those of the diagonal terms |diha| and |aiha|. In Figure 2, we also show k(t) calculated using the electronic acceptor population operator pa = |aiha| instead of the acceptor population operator h defined in Eq. (3) to calculate the product state population. It can be seen that the calculated rate constants agree very well, although the initial relaxations of k(t) using the two definitions are very different. This alternative definition of the product state population using the pa operator also works well for each components of the initial density matrix (results not shown), which indicates that after the initial relaxation, both pa and h can be used to describe the product population. Calculating ET rate constants using the acceptor state population operator pa has another advantage that the reactive flux method can also be applied to ET reaction in the inverted region. Figure 3 shows the contribution to k(t) from the diagonal |dihd| and |aiha| terms with ˙ It can be seen that, for the |dihd| term, positive and negative initial collective velocity X. transmission coefficients from the trajectories with positive and negative initial velocities are both small, while the contribution from those with positive initial velocity is slightly larger. For the |aiha| term, the small transmission coefficient shown in Figure 2 is actually due to the cancellation of contributions from trajectories with positive and negative initial collective velocities. Such behavior agrees with previous discussions in Section 2.2. We then diagonalize the mixed quantum-classical Hamiltonian HS + X ∗ σz at the transition state and calculate contributions to k(t) from the adiabatic states |1i and |2i. The results are shown in Figure 4. We can see that contribution from the upper state |2i is negative, which means that if the electronic state is initially on the adiabatic state |2i and the

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

initial collective velocity is positive, the final state tends to be in the donor state comparing to the case with negative initial collective velocity. Since contributions from the diagonal terms in the diabatic basis are all positive, this behavior is due to contributions from the off-diagonal |diha| and |aihd| terms. We note that the contributions of off-diagonal terms shown in Figure 2 are weighted by the equilibrium values of the initial density matrix ρ0da and ρ0ad , which is proportional to V in the non-adiabatic limit. To further investigate the effects of the off-diagonal terms, we calculate the unweighted contributions to k from the off-diagonal density operator |diha| with different strength of electronic coupling Vd (k(t) calculated from the unweighted |diha| term is actually negative, so the absolute values are used), and the results are shown in Figure 5, where they are also compared with the total ET rate constants. By investigating the different slopes of two curves in the log-log plot, it can be seen that the transmission coefficient contributed from the unweighted off-diagonal term is actually proportional to the coupling V in the non-adiabatic limit, in contrast to the V 2 scaling of the total rate constant. In the remaining part of this section, we use the MQCL results to test the validity of the Ehrenfest dynamics in calculating the ET rate constants using the reactive flux formalism. In a previous paper, 16 we have shown that, when starting from the equilibrium nuclear configuration on the donor state, the Ehrenfest dynamics is capable to give correct ET rate constants in a wide range of parameter regimes. However, in the reactive flux method, the initial density matrix include both diagonal and off-diagonal terms, and it is not clear whether the mean field description in the Ehrenfest dynamics is sufficient to obtain the correct ET rate constants. Figure 6 shows the result of total k(t), where the initial density matrix in Eq. (10) is used in the calculation with the Ehrenfest dynamics. Results using both the h and pa operators to define the product population are presented. It can be seen that, unlike the case of MQCL simulation, results from the two different definition no longer agree with each other in the Ehrenfest dynamics, and neither of the two definitions give good agreement with the exact

12


Page 12 of 30

Page 13 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


MQCL result. We can also decompose the initial density matrix into the diagonal and off-diagonal terms, and calculate their contributions to the total ET rate. Unlike the MQCL method, the Ehrenfest dynamics is not additive with regard to the initial density matrix. Figure 7 shows the results for contributions from the individual diagonal (only the |dihd| term is shown), and off-diagonal (|diha| plus |aihd|) terms. We can see that, when using the Heaviside function h to define the product population, the Ehrenfest dynamics still fails badly to reproduce the MQCL result. However, when using the pa operator to define the product state population, results from the Ehrenfest dynamics agree reasonably well with the MQCL results. Figure 8 shows the ET rate calculated using the Ehresfest dynamics in combination with the reactive flux method, by summing up the results using the separated |dihd|, |aiha|, and |diha| + |aihd| terms of the density matrix, and using the pa operator to define the product population. It can be seen that the approximate Ehrenfest dynamics give reasonably good results when comparing to the numerically exact MQCL method.

4

Conclusions

The reactive flux formalism for ET (and non-adiabatic PT) reaction rate constants have been derived and used in many previous works. 20,21,23–26,31,50 However, exactly following the non-adiabatic dynamics to calculate the flux-population correlation function is still not an easy task. Many previous works have employed transmission coefficients based on the Landau-Zener formula 8,9 and the Zhu-Nakamura theory, 10,11 as well as mixed quantumclassical dynamics such as the surface hopping method 20,21,23 and the semiclassical method based on the mapping Hamiltonian. 31 In this work, we have calculated ET rate constants using the MQCL method which is exact for the spin-Boson model. Together with previous studies by Kapral and coworkers using the MQCL method in the adiabatic basis, 24–26 it is shown that the MQCL method can be combined effectively with the reactive flux formalism

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to calculate non-adiabatic ET rate constants. Contributions from the different terms of the initial density matrix at the transition state using the MQCL equation in the diabatic basis are investigated. It is found that the off-diagonal terms of the initial density matrix play an important role in determining the transmission coefficient. This is more clearly demonstrated when calculating the unweighted contributions from the off-diagonal terms to the rate constant, which is found to be proportional to the electronic coupling Vd in the non-adiabatic limit, and results in negative contribution to the total rate constant from the upper adiabatic state. Due to the limit of current numerical methods to propagate the MQCL equation, most simulations performed in this work is in the non-adiabatic regime with Vd ≪ kB T . In future studies, it would also be interesting to investigate the effect of off-diagonal terms beyond the non-adiabatic limit, where Vd is larger or the temperature is lower, and the ground adiabatic state plays a more important role in determining the total transmission coefficient. We have also compared the results with those from the Ehrenfest dynamics using the reactive flux method. The Ehrefest dynamics is found to be incorrect when using the total mixed density matrix as the initial state for the electronic degree of freedom. However, when computing contributions from the individual elements of the initial density matrix and using the electronic acceptor population operator pa to define the product population, the Ehrenfest dynamics may also give reasonably accurate results. This also agrees with the previous finding that, direct calculation of the ET rate using Ehrenfest dynamics starting from equilibrated nuclear configurations on the donor state can result in correct rate constants. 16 Currently, simulations using the MQCL equation still suffer from the problem of numerical convergence at longer times, yet it is still possible to develop new approximate methods based on the MQCL equation to calculate non-adiabatic reaction rate constants, which could be an interesting topic in future studies.

14


Page 14 of 30

Page 15 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgement This work is supported by NSFC (Grant No. 21290194), the 973 program (Grant No. 2013CB933501), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB12020300).

A

Calculation of the transition state probability

The Z ∗ /Z term defined in Eq. (9) can be written as

Z = Z ∗

*

δ

X

cj xj − X ∗

j

!+

=

h P i ∗ Tr e−βH δ c x − X j j j Tre−βH

,

(14)

where X ∗ is the dividing surface between the reactant and product states. By defining "

−βH

Z(X) = Tr e

δ

X

cj xj − X

j

!#

,

(15)

Z ∗ /Z in Eq. (14) can be calculated as Z(X ∗ ) Z∗ = R +∞ . Z dXZ(X)

(16)

−∞

Detailed calculation of Z(X) is very similar to that in a previous study of PT in a double well model coupled to a harmonic bath. 50 So only the final result is given below. It is found that

where

X2 TrS e−β(Hs +Xσz ) , Z(X) = C1 exp − λkB T λ=

X 2c2j mj ωj2 j

(17)

(18)

is the reorganization energy. The Z ∗ /Z term can thus be calculated using Equations (16-18).

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 30

To calculate the free energy surface shown in Fig. 2, we further define X2 Z1,2 (X) = C1 exp − e−βE1,2 , λkB T

(19)

where the adiabatic state energies E1 and E2 are the eigenvalues of the mixed quantumclassical Hamiltonian Hs + Xσz . The corresponding free energy surfaces are then defined as 1 W1,2 (X) = − ln Z1,2 (X) . β

B

(20)

Sampling the initial nuclear positions and velocities

For the mixed quantum-classical Boltzmann distribution in Eq. (10), we can see that the first term is just an equilibrated Boltzmann density matrix for the electronic degrees of freedom, and the second term is the Boltzmann operator of the bare bath (e−βHB ) times the delta function δ(X − X ∗ ), which can be written as −β

e

P

j

pj 2 2mj

+ 21 mj ωj2 x2j

δ

X

cj xj − X ∗

j

!

.

(21)

The initial velocities are thus sampled directly from the classical Boltzmann distribution. To sample the initial nuclear positions x, we first use the following relation (assuming that we have N bath modes) X

cj xj = X ∗

(22)

j

to write xN as a function of {x1 , · · · , xN −1 }. The distribution of the bath positions in Eq. (21) is thus a Gaussian function of N − 1 variables. We can thus sample the N − 1 position variables {x1 , · · · , xN −1 } firstly, and obtain xN through Eq. (22).

16


Page 17 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References (1) May, V.; K¨ uhn, O. Charge and Energy Transfer Dynamics in Molecular Systems, 3rd ed.; Wiley-VCH: Weinheim, 2011. (2) Domcke, W.; Yarkony, D.; K¨oppel, H. Conical Instersections: Electronic Structure, Dynamics and Spectroscopy; World Scientific: Singapore, 2004. (3) Marcus, R. A. On the Theory of Oxidation-Reduction Involving Electron Transfer. J. Chem. Phys. 1956, 24, 966–978. (4) Marcus, R. A. Electron Transfer Reactions in Chemistry. Theory and Experiment. Rev. Mod. Phys. 1993, 65, 599–610. (5) Weiss, S. Fluorescence Spectroscopy of Single Biomolecules. Science 1999, 283, 1676. (6) Landau, L. On the Theory of Transfer of Energy at Collisions II. Physik. Z. Sow. 1932, 2, 46. (7) Zener, C. Non-adiabatic Crossing of Energy Levels. Proc. R. Soc. London Ser. A 1932, 137, 696–702. (8) Cline, R. E.; Wolynes, P. G. Stochastic Dynamic Models of Curve Crossing Phenomena in Condensed Phases. J. Chem. Phys. 1987, 86, 3836. (9) Straub, J. E.; Berne, B. J. A Statistical Theory for the Effect of Nonadiabatic Transitions on Activated Processes. J. Chem. Phys. 1987, 87, 6111. (10) Zhao, Y.;

Milnikov, G.;

Nakamura, H. Evaluation of Canonical and Micro-

canonical Non-adiabatic Raction Rate Constants by Using Zhu-Nakamura Formulas. J. Chem. Phys. 2004, 121, 8854–8860. (11) Zhao, Y.; Liang, W. Z.; Nakamura, H. Semi-classical Treatment of Thermally Activated

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 30

Electron Transfer in the Intermediate to Strong Electronic Coupling Regime Under the Fast Dielectric Relaxation. J. Phys. Chem. A 2006, 110, 8204–8212. (12) Zhu, C.; Nakamura, H. Theory of Nonadiabatic Transition for General Two-state Curve Crossing Problems. I. Nonadiabatic Tunneling Case. J. Chem. Phys. 1994, 101, 10630. (13) Landry, B. R.; Subotnik, J. E. Communication: Standard Surface Hopping Predicts Incorrect Scaling for Marcus’ Golden-rule Rate: The Decoherence Problem Cannot be Ignored. J. Chem. Phys. 2011, 135, 191101–191104. (14) Landry, B. R.; Subotnik, J. E. How to Recover Marcus Theory With Fewest Switches Surface Hopping: Add Just a Touch of Decoherence. J. Chem. Phys. 2012, 137, 22A513. (15) Landry, B. R.; Falk, M. J.; Subotnik, J. E. Communication: The Correct Interpretation of Surface Hopping Trajectories: How to Calculate Electronic Properties. J. Chem. Phys. 2013, 139, 211101. (16) Xie, W.; Bai, S.; Zhu, L.; Shi, Q. Calculation of Electron Transfer Rates Using Mixed Quantum Classical Approaches: Nonadiabatic Limit and Beyond. J. Phys. Chem. A 2013, 117, 6196–6204. (17) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (18) H¨anggi, P.; Talkner, P.; Borkovec, M. Reaction Rate Theory: Fifty Years after Kramers. Rev. Mod. Phys. 1990, 62, 251. (19) Nitzan, A. Chemical Dynamics in Condensed Phases; Oxford University Press: New York, 2006. (20) Hammes-Schiffer, S.; Tully, J. C. Nonadiabatic Transition State Theory and Multiple Potential Energy Surface Molecular Dynamics of Infrequent Events. J. Chem. Phys. 1995, 103, 8528–8537. 18


Page 19 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(21) Kim, S. Y.; Hammes-Schiffer, S. Hybrid Quantum/classical Molecular Dynamics for a Proton Transfer Reaction Coupled to a Dissipative Bath. J. Chem. Phys. 2006, 124, 244102. (22) Topaler, M.; Makri, N. Quamtum Rates for a Double Well Coupled to a Dissipative Bath: Accurate Path Integral Results and Comparison with Approximate Theories. J. Chem. Phys. 1994, 101, 7500. (23) Jain, A.; Subotnik, J. E. Surface Hopping, Transition State Theory, and Decoherence. II. Thermal Rate Constants and Detailed Balance. J. Chem. Phys. 2015, 143, 134107. (24) Sergi, A.; Kapral, R. Quantum-classical Dynamics of Nonadiabatic Chemical Reactions. J. Chem. Phys. 2003, 118, 8566. (25) Sergi, A.; Kapral, R. Nonadiabatic Reaction Rates for Dissipative Quantum-Classical Systems. J. Chem. Phys. 2003, 119, 12276. (26) Kim, H.; Kapral, R. Nonadiabatic Quantum-Classical Reaction Rates with Quantum Equilibrium Structure. J. Chem. Phys. 2005, 123, 194108–194117. (27) Hanna, G.; Kapral, R. Quantum-classical Liouville Dynamics of Nonadiabatic Proton Transfer. J. Chem. Phys. 2005, 122, 194108. (28) Kim, H.; Hanna, G.; Kapral, R. Analysis of Kinetic Isotope Effects for Nonadiabatic Reactions. J. Chem. Phys. 2006, 125, 084509. (29) Kim, H.; Kapral, R. Solvation and Proton Transfer in Polar Molecule Nanoclusters. J. Chem. Phys. 2006, 125, 234309. (30) 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.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 30

(31) Ka, B. J.; Thompson, W. H. Nonadiabatic Effects on Proton Transfer Rate Constants in a Nanoconfined Solvent. J. Phys. Chem. B 2010, 114, 7535. (32) Azzouz, H.; Borgis, D. A Quantum Molecular-dynamics Study of Proton-transfer Reactions along Asymmetrical H Bonds in Solution. J. Chem. Phys. 1993, 98, 7361. (33) Miller, W. H. Electronically Nonadiabatic Dynamics Via Semiclassical Initial Value Methods. J. Phys. Chem. A 2009, 113, 1405–1415. (34) Nan, G.-J.; Zheng, R.-H.; Shi, Q. Mixed Quantum-classical Approaches to Calculating Charge Transfer Rate Constants: Applications to Realistic Systems. Acta Phys. Chim. Sinica 2010, 26, 1755. (35) Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M.; Garg, A.; Zwerger, W. Dynamics of the Dissipative Two-level System. Rev. Mod. Phys. 1987, 59, 1. (36) Weiss, U. Quantum Dissipative Systems, 2nd ed.; World Scientific: London, 1999. (37) Yamamoto, T. Quantum Statistical Mechanical Theory of the Rate of Exchange Chemical Reactions in the Gas Phase. J. Chem. Phys. 1960, 33, 281. (38) Gerasimenko, V. I. Dynamical Equations of Quantum-Classical: Systems. Theor. Math. Phys. 1982, 50, 49–55. (39) Kapral, R.; Ciccotti, G. Mixed Quantum-Classical Dynamics. J. Chem. Phys. 1999, 110, 8919–8929. (40) Kapral, R. Quantum-Classical Dynamics in a Classical Bath. J. Phys. Chem. A 2001, 105, 2885–2889. (41) Ando, K.; Santer, M. Mixed Quantum-Classical Liouville Molecular Dynamics without Momentum Jump. J. Chem. Phys. 2003, 118, 10399–10406.

20


Page 21 of 30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(42) Shi, Q.; Geva, E. A Derivation of the Mixed Quantum-Classical Liouville Equation from the Influence Functional Formalism. J. Chem. Phys. 2004, 121, 3393–3404. (43) Tanimura, Y.; Mukamel, S. Multistate Quantum Fokker-planck Approach to Nonadiabatic Wave-packet Dynamics in Pump probe Spectroscopy. J. Chem. Phys. 1994, 101, 3049–3061. (44) Horenko, I.; Salzmann, C.; Schmidt, B.; Sch¨ utte, C. Quantum-Classical Liouville Approach to Molecular Dynamics: Surface Hopping Gaussian Phase-Space Packets. J. Chem. Phys. 2002, 117, 11075–11088. (45) Santer, M.; Manthe, U.; Stock, G. Quantum-classical Liouville Description of Multidimensional Nonadiabatic Molecular Dynamics. J. Chem. Phys. 2001, 114, 2001–2012. (46) Mac Kernan, D.; Ciccotti, G.; R.Kapral, Surface-Hopping Dynamics of a Spin-Boson System. J. Chem. Phys. 2002, 116, 2346–2353. (47) Kapral, R. Progress in The Theory of Mixed Quantum-Classical Dynamics. Annu. Rev. Phys. Chem. 2006, 57, 129–157. (48) Shi, Q.; Geva, E. A Self-consistent Treatment of Electron Transfer in the Limit of Strong Friction Via the Mixed Quantum Classical Liouville Method. J. Chem. Phys. 2009, 131, 034511. (49) Mac Kernan, D.; Ciccotti, G.; R.Kapral, Trotter-Based Simulation of QuantumClassical Dynamics. J. Phys. Chem. B 2008, 112, 424–432. (50) Xie, W.; Xu, Y.; ; Zhu, L.; Shi, Q. Mixed Quantum Classical Calculation of Proton Transfer Reaction Rates: From Deep Tunneling to over the Barrier Regimes. J. Chem. Phys. 2014, 140, 174105.

21



0.06

Free Energy (a.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 30

0.004 0.002

0.04

0 -0.002

0

0.002

0.02

0

-0.04

-0.02

0

0.02

0.04

X (a.u)

Figure 1: Free energy profiles as a function of the collective reaction coordinate X, for the spin-Boson model studied in Section 3. The solid black line is for the adiabatic ground state, the dashed red line is for the adiabatic excited state.

22


Page 23 of 30

MQCL, Total MQCL, |d>

Mixed Quantum-Classical Study of Nonadiabatic Curve Crossing in

Recommend Documents