Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (1): 12-26

#### The article information

Liang-hui Gao, Bin-bin Xie, Wei-hai Fang

Theories and Applications of Mixed Quantum-Classical Non-adiabatic Dynamics

Chinese Journal of Chemical Physics, 2018, 31(1): 12-26

http://dx.doi.org/10.1063/1674-0068/31/cjcp1712234

### Article history

Received on: December 6, 2017
Accepted on: January 6, 2018
Theories and Applications of Mixed Quantum-Classical Non-adiabatic Dynamics
Liang-hui Gao, Bin-bin Xie, Wei-hai Fang
Dated: Received on December 6, 2017; Accepted on January 6, 2018
Key Laboratory of Theoretical and Computational Photochemistry, Ministry of Education, College of Chemistry, Beijing Normal University, Beijing 100875, China
*Author to whom correspondence should be addressed. Liang-hui Gao, E-mail:lhgao@bnu.edu.cn; Wei-hai Fang, E-mail:fangwh@bnu.edu.cn
Abstract: Electronically non-adiabatic processes are essential parts of photochemical process, collisions of excited species, electron transfer processes, and quantum information processing. Various non-adiabatic dynamics methods and their numerical implementation have been developed in the last decades. This review summarizes the most significant development of mixed quantum-classical methods and their applications which mainly include the Liouville equation, Ehrenfest mean-field, trajectory surface hopping, and multiple spawning methods. The recently developed quantum trajectory mean-field method that accounts for the decoherence corrections in a parameter-free fashion is discussed in more detail.
Key words: Non-adiabatic dynamics    Mean-field    Surface hopping    Decoherence
Ⅰ. INTRODUCTION

The motion and energy and momentum exchanges of a chemical system are governed by the laws of quantum mechanics (QM). It is essential to directly solve the Schrödinger equation of both the electrons and nuclei. Nevertheless, because the mass of electron is lighter than nucleus by 3 orders, it is natural to introduce separation between time/energy scales of electrons and nuclei, the so-called Born-Oppenheimer (BO) approximation [1], which distinguishes the rapidly moving electrons from the slowly moving nuclei. The BO representation leads to the formation of electronic adiabatic eigenstates with fixed nuclear position. The nuclear motion is then governed by a single adiabatic BO potential energy surface (PES) [2], which describes the variation of the electronic energy with changes in the nuclear geometry. The full-dynamics simulations can be carried on-the-fly by assuming nuclear motion with forces obtained from the gradient of the PES using ab initio quantum chemistry.

The BO approximation is appropriate for many chemical reactions occurring entirely on a single electronic state and not involving proton or electron transfer. However, there are a great number of chemical events that cannot be described within the BO approximation [3-6]. For example, in many photochemical processes, the nuclear dynamics cannot be adequately modeled on a single PES because for some nuclear configurations the separation between PESs becomes comparable to the nuclear energy scale or even vanishes (the case referred as conical intersections). As a result, the electronic states change during the dynamic events. By energy conservation, every electronic transition downward must be associated with an increase in the nuclear thermal energy. This energy can in turn push the molecule to move toward high energy geometries, where new electronic state crossing might be found, and the cycle continues. The electronically non-adiabatic processes are also essential parts of collisions of electronically excited species, chemiluminescent reactions, recombination reactions, heterolytic dissociations, and electron transfer processes [6-11]. Non-adiabatic dynamics also play intimate roles in quantum information processing with atomic, molecular, or solid-state qubits [12]. Proper incorporation of non-adiabatic electronic transition in theoretical and computational methods is thus crucial for understanding and predicting the outcomes of these classes of important chemical processes.

A standard picture of non-adiabatic process was introduced in early 1930's when quantum mechanism theory was just established [13, 14]. It involves transition between diabatic surfaces corresponding to coupled reactant and product sites. The Hamiltonian for such coupled two-states system [15] is

 $\hat H = \frac{{{\bf{P}}_s^2}}{{2m}}I + {\rm{ }}\left( {\frac{{{V_1}(s)}}{{{V_d}(s)}}\frac{{{V_d}(s)}}{{{V_2}(s)}}} \right)$ (1)

Here $I$ is the unit matrix and $s$ is a nuclear reaction coordinate. The diabatic potential matrix is a full matrix in which the diagonal elements are the potential energies for the elastic motion in two particular states, while the off-diagonal elements provide the coupling between these two states. A sketch of the diabatic potential surfaces is given in FIG. 1; the two curves are assumed to cross at $s$=0. By diagonalizing the potential energy part of $\hat H$, two adiabatic energy surfaces are obtained,

 $V_ \pm \left( s \right) = \frac{1}{2}\left( {{V_1}\left( s \right) + {V_2}\left( s \right)} \right) \pm \\ ~~~~~~~~~~{}\frac{1}{2}\sqrt {{{\left( {{V_1}\left( s \right) - {V_2}\left( s \right)} \right)}^2} + 4V_d^2\left( s \right)}$ (2)
 FIG. 1 Schematic diagram of the diabatic (black solid lines) and adiabatic (red dashed lines) potential surfaces.

The system exhibits an "avoid crossing" between the two adiabatic surfaces, see FIG. 1. If the energy gap between the PESs is sufficiently large, then adiabatic description of the reaction under BO approximation is adequate. The deviation from BO single-surface dynamics occurs in the region of configuration spaces where the two PESs nearly cross. Under the assumption that near the conical crossing point ($s$=0) the velocity $v$ of the particle is constant, the two crossing diabatic potential curves are linear, and the off-diagonal coupling $V_d$ is a constant. Landau and Zener (LZ) derived the formula for the probability for the transition passage from one state to another [13, 14],

 ${P_{{\rm{LZ}}}} = \exp \left( {\frac{{ - 2\pi V_d^2}}{{\hbar v\left| {V_1'\left( 0 \right) - V_2'\left( 0 \right)} \right|}}} \right)$ (3)

Here, $V_1'$(0) and $V_2'$(0) are the first derivatives of the diagonal potentials. The exponent in Eq.(3), called Massey parameter, determines whether the transition is adiabatic (if it is much greater than unity) or completely non-adiabatic (if it is much less than unity).

The LZ-formula provides a good approximation for non-adiabaticity only in the vicinity of the crossing point. The missing of the nuclear kinetic energy in LZ model obstacles the first principles description of chemical reactivity. Therefore, the development of ab initio molecular dynamics that can properly involve the non-adiabatic processes must meet the following requirements: first, the electronic problem must be solved efficiently and accurately for both the ground and excited states; second, the quantum mechanical character of the nuclear dynamics must be addressed because at least two electronic state will be involved during the dynamics.

Much progress of non-adiabatic molecular dynamics has been made since 1970's [4, 5], including full quantum mechanical methods where both electrons and nuclei are described via wave packets, mixed quantum-classical (MQC) approaches where the electrons are treated as quantum particle while the nuclei are treated as classical particles, and semi-classical method that use a semi-classical description for all the degree of freedom of electrons and nuclei. Time-dependent Hartree and multi-configuration time-dependent Hartree are fully quantum mechanical methods that propagate wave packets for nuclei [16-25]. However, a full quantum solution of the Schrödinger equation is only feasible for systems with a few degrees of freedom (DOF) and the Hamiltonians have special forms. It is not practical for treating non-adiabatic dynamics in real system. For larger systems with thousands of DOF, it is more common to separate the system to a quantum subsystem coupled to a classical subsystem. In this situation, the light quantum particles are represented by quantum-mechanical description, while the heavy classical particles are represented by a trajectory description. Mixed quantum-classical Liouville (QCL) equation in the form of density matrix is a promising method [26-31]. A variety of quantum-classical trajectory approaches have also been proposed, including Ehrenfest mean-field (EMF) method [32-38], trajectory surface hopping (TSH) method [39-68], quantum-trajectory mean-field (QTMF) method [69], and multiple spawning (MS) method [70-75]. In this review, we summarize the most significant development of the mixed quantum-classical methods and their applications. Because the coupling of quantum and classical DOF in mixed quantum-classical approaches is not completely consistent, an alternative strategy is to use a semi-classical description for all the degrees of freedom [76-81]. As proposed by Meyer and Miller, the electronic state is characterized by a pair of classic harmonic oscillator action-angle variable, which can be transformed into Cartesian coordinates and momenta for the actual numerical trajectory calculation. We suggest the readers to refer to Refs.[76-81] for a detailed discussion and review of semi-classical method.

Ⅱ. NON-ADIABATIC DYNAMICS THEORIES AND APPLICATIONS A. Mixed quantum-classical Liouville equation

For condensed phase system or vibration motions of molecules in liquid with thousands of degrees of freedom., it is not feasible to attempt a full quantum solution of the Schrödinger equation. Consequently, one is led to considering dynamics of quantum subsystem coupled to a classical subsystem. In this situation, the light quantum particles are represented by a quantum mechanical description, while the heavy classical particles are represented by a trajectory description [26, 27]. Suppose that the quantum-classical system is composed of an $n$-particle QM subsystem and an $N$-particle classical subsystem with coordinate operator $\bf\hat q$ and $\bf\hat Q$, momentum operator $\bf\hat p$ and $\bf\hat P$, and mass $m$ and $M$, respectively, then the Hamiltonian is

 $\hat H=\frac{{{{\bf{\hat P}}^2}}}{{2M}} + \frac{{{{\bf{\hat p}}^2}}}{{2m}} + \hat V\left( \bf{\hat q, \hat Q} \right)$ (4)

The evolution of the density matrix of the full system is given by quantum Liouville equation

 $\frac{{\partial \hat \rho }}{{\partial t}} = - \frac{i}{\hbar }\left[{\hat H, \hat \rho } \right]$ (5)

To treat the nuclear degree of freedom as classic, a partial Wigner transformation with respect to the Q-coordinates is essentially performed on the density matrix $\hat\rho(t)$,

 ${\hat \rho _W}\left( {\bf{R}}, {\bf{P}} \right) = \frac{1}{\left( 2\pi \hbar\right)^{3N}}\int{\rm{d}}{\bf{z}}\textrm{e}^{i\bf{P\cdot z}/\hbar }\\ ~~~~~~~~~~~~~~~~~~~~{}\bigg\langle{\bf{R}} - \frac{\bf{z}}{2}\left|\hat \rho \right|{\bf{R}} + \frac{\bf{z}}{2} \bigg\rangle$ (6)

In the limit of $M$$\gg$$m$, Kapral and Ciccotti (KC) [27] obtained the evolution equation of $\hat\rho_W({\bf{R}}, {\bf{P}}, t)$,

 $\frac{\partial\hat\rho_W({{\bf{R}}}, {{\bf{P}}}, t)}{\partial t} = - \frac{i}{\hbar}\left[\hat H_W({\bf{R}}, {\bf{P}}), \hat\rho_W({\bf{R}}, {\bf{P}}, t) \right] +\\~~~~~~~~~~~~~~~~~~~~ {}\frac{1}{2}\Big(\left\{\hat H_W({\bf{R}}, {\bf{P}}), \hat\rho_W({\bf{R}}, {\bf{P}}, t)\right\}-\\~~~~~~~~~~~~~~~~~~~~ {}\left\{\hat\rho_W({\bf{R}}, {\bf{P}}, t), \hat H_W({\bf{R}}, {\bf{P}})\right\}\Big)$ (7)

Here

 ${\hat H_W}\left({\bf{R}}, {\bf{P}} \right) = \frac{{{{\bf{P}}^2}}}{{2M}} + \frac{{{{\bf\hat p}^2}}}{{2m}} + {\hat V_W}\left( {\bf\hat q, {\bf{R}}} \right)$ (8)

Compared to Eq.(4), the remove of hat from $\bf\hat P$ and $\bf\hat Q$ implies that the $M$-particles are treated classically. The Poisson bracket is defined by

 $\left\{ {\hat A\left({\bf{R}}, {\bf{P}} \right), \hat B\left({\bf{R}}, {\bf{P}}\right)} \right\} = \frac{{\partial \hat A}}{{\partial{\bf{R}}}} \cdot \frac{{\partial \hat B}}{{\partial{\bf{P}}}} - \frac{{\partial \hat A}}{{\partial {\bf{P}}}} \cdot \frac{{\partial \hat B}}{{\partial {\bf{R}}}}$ (9)

Eq.(7) is referred as the quantum-classical Liouville equation.

To solve the evolution equation of $\hat\rho_W({\bf{R}}, {\bf{P}}, t)$, a convenient basis of adiabatic states $|\alpha;{\bf{R}}\rangle$ is selected, which is defined by

 ${\hat h_W}\left( {\bf{R}} \right)|\alpha ; {\bf{R}}\rangle = {E_\alpha }\left( {\bf{R}} \right)|\alpha ;{\bf{R}}\rangle$ (10)

The Hamiltonian $\hat h_W({\bf{R}})$ is for the quantum subsystem with fixed values of the classic coordinates ${\bf{R}}$, which has form

 ${\hat h_W}\left(\bf R \right) = \frac{{{{\bf\hat p}^2}}}{{2m}} + {\hat V_W}\left( {\bf\hat q, {\bf{R}}} \right)$ (11)

Taking the matrix elements of $\hat\rho_W$ as

 $\rho _W^{\alpha \alpha '}\left( {\bf{R}}, {\bf{P}} \right) = \langle\alpha ;{\bf{R}}|{\hat \rho _W}\left( {\bf{R}}, {\bf{P}} \right)|\alpha ';\bf R\rangle$ (12)

the adiabatic representation of the equation of motion for the QC system is

 $\frac{{\partial \rho _W^{\alpha \alpha '}\left( {\bf{R}}, {\bf{P}} \right)}}{{\partial t}} = - i{\omega _{\alpha \alpha '}}\left( {\bf{R}} \right)\rho _W^{\alpha \alpha '}\left({\bf{R}}, {\bf{P}} \right) +\\~~~~~~~~~~~~ {}\sum\limits_{\alpha'' } \frac{\bf P}{M} \cdot \left( {\rho _W^{\alpha \alpha ''}{\bf{d}_{\alpha'' \alpha '}} - {\bf{d}_{\alpha \alpha'' }}\rho _W^{\alpha'' \alpha '}} \right)-\\~~~~~~~~~~~~ \frac{{\bf{P}}}{M} \cdot \frac{\partial \rho _W^{\alpha \alpha '}}{\partial {\bf{R}}} - \frac{1}{2}\sum\limits_{\alpha''} \Big({\bf F}_W^{\alpha \alpha ''} \cdot \frac{\partial \rho _W^{\alpha'' \alpha '}}{\partial{\bf{P}}} +\\~~~~~~~~~~~~ {}\frac{\partial \rho _W^{\alpha \alpha'' }}{\partial{\bf{P}}}\cdot F_W^{\alpha'' \alpha '} \Big)$ (13)

Here,

 ${\omega _{\alpha \alpha '}}\left({\bf{R}} \right) = \frac{1}{\hbar }\left( {{E_\alpha }\left( {\bf{R}} \right) - {E_{\alpha '}}\left( {\bf{R}} \right)} \right)$ (14)

is the transition frequency between two electronic states,

 $\bf{d}_{\alpha \alpha '} = \left\langle\alpha ;{\bf{R}} \left|\frac{\partial }{\partial {\bf{R}} }\right|\alpha ';{\bf{R}} \right\rangle$ (15)

is the non-adiabatic coupling matrix element, and

 ${\bf F}_W^{\alpha \alpha '} = - \left\langle\alpha ;{\bf{R}}\left|\frac{\partial\hat{V}_W\left({\bf\hat{q}, R} \right)}{\partial{\bf{R}}}\right|\alpha ';{\bf R}\right\rangle$ (16)

is the Hellmann-Feynman force.

Ando presented an alternative QCL formulation by taking a different route, namely, by first taking the matrix elements over the adiabatic electronic basis, and then the partial Wigner transformation for the nuclear coordinates [29], he obtained

 $\frac{{\partial \rho _W^{\alpha \alpha '}\left({\bf{R}}, {\bf{P}} \right)}}{{\partial t}} = - i{\omega _{\alpha \alpha '}}\left( {\bf{R}} \right)\rho _W^{\alpha \alpha '}\left( {\bf{R}}, {\bf{P}} \right) +\\~~~~~~~~~~~~ {} \sum\limits_{\alpha''}\frac{{\bf{P}}}{M}\cdot\left( {\rho _W^{\alpha \alpha ''}{\bf{d}_{\alpha'' \alpha '}} - {\bf{d}_{\alpha \alpha'' }}\rho _W^{\alpha'' \alpha '}} \right)- \\~~~~~~~~~~~~ {}\frac{{\bf{P}}}{M} \cdot\frac{{\partial \rho _W^{\alpha \alpha '}}}{{\partial {\bf{R}}}} -\frac{1}{2}\left( {{{\bf F}_\alpha }\left( {\bf{R}} \right) + {{\bf F}_\beta }\left( {\bf{R}} \right)} \right) \cdot\\~~~~~~~~~~~~ {} \frac{{\partial \rho _W^{\alpha \alpha '}\left( {\bf{R}}, {\bf{P}} \right)}}{{\partial {\bf{P}}}} +{\cal O}\left( \hbar \right)$ (17)

Here $\bf{F}_\alpha({\bf{R}})$=$- \partial{E_\alpha}({{\bf{R}}})/\partial {{\bf{R}}}$ is the adiabatic force. Compared to the KC equation, the Ando-QCL equation is in absence of the off-diagonal Hellmann-Feynman forces.

The QCL equation can be solved numerically by employing stochastic algorithm using a Monte Carlo sampling of local classical trajectories [28, 30]. The formal solution of the QCL equation is

 $\rho \left( t \right) = {{\rm{e}}^{ - i\left( {{{\cal L}_0} + {{\cal L}^{\textrm{QC}}}} \right)t/\hbar }}\rho \left( 0 \right)$ (18)

where the Liouville operator ${{\cal L}_0}$=${{\cal L}^\textrm{Q}}$+${{\cal L}^\textrm{C}}$ accounts for the pure quantum-mechanical and pure classical time evolution of the density matrix with

 ${\cal L}_{\alpha \alpha ', \beta \beta '}^\textrm{Q} = - \frac{i}{\hbar }\left( {{E_\alpha }\left( \bf R \right) - {E_{\alpha '}}\left(\bf R \right)} \right){\delta _{\alpha \beta }}{\delta _{\alpha '\beta '}}$ (19)
 ${\cal L}_{\alpha \alpha ', \beta \beta '}^\textrm{C} =\bigg[{\frac{\bf P}{M} \cdot \frac{\partial }{{\partial\bf R}}-\frac{1}{2}\left( {{{\bf F}^{\alpha \alpha }}\left(\bf R \right) + {{\bf{F}}^{\beta \beta }}\left(\bf R \right)} \right)}\cdot\\~~~~~~~~~~~~~~~~~ {}{\frac{\partial }{{\partial\bf P}}}\bigg]{\delta _{\alpha \beta }}{\delta _{\alpha '\beta '}}$ (20)

The QC part of the Liouville operator reads

 ${\cal L}_{\alpha \alpha ', \beta \beta '}^{\textrm{QC}} = - \frac{\bf P}{M} \cdot {{\bf d}_{\alpha \beta }}\left( {1 + \frac{1}{2}{{\bf S}_{\alpha \beta }} \cdot \frac{\partial }{{\partial\bf P}}} \right){\delta _{\alpha '\beta '}}\cdot\\~~~~~~~~~~~~~~~~~ {}\left( {1 - {\delta _{\alpha \beta }}{\delta _{\alpha '\beta '}}} \right) + {\rm{h}}.{\rm{c}}.$ (21)
 ${S_{\alpha \beta }} = \frac{{{E_\alpha }\left(\bf R \right) - {E_\beta }\left(\bf R \right)}}{{\displaystyle\frac{\bf P}{M} \cdot {\bf d_{\alpha \beta }}}}{{\bf d}_{\alpha \beta }}$ (22)

By using the first-order Trotter approximation

 $\rho \left( {t + \delta t} \right) = {{\rm{e}}^{ - \frac{i}{\hbar }{{\cal L}_0}\delta t}}\left( {1 - \frac{i}{\hbar }{{\cal L}^{\textrm{QC}}}\delta t} \right)\rho (t) + {\cal O}(\delta {t^2})$ (23)

The equation of motion for a small time step $\delta t$ can be solved in two processes: (ⅰ) classical propagation and phase accumulation on the PES described by ${\cal L}_0$; (ⅱ) transition between coupled PES induced by ${\cal L}^{\textrm{QC}}$. Then a phase-space integral is performed to evaluate the time-dependent observables. In the KC-QCL equation, the off-diagonal term ($S_{\alpha\beta}$) has been attributed to a meaning of "momentum-jump" associated with the non-adiabatic transition in the trajectory implementation. However, numerical calculation showed that the momentum-jump operation resulted in instability around the classical turning points of the nuclear motion [28, 31].

The QCL approach provides a correct quantum-classical description of the electronic coherence. The limitation of this approach is that a huge number of trajectories in the phase space are needed. The sampling effort increases exponentially in time. It also has an implementation problem of numerical instability [28]. Therefore, this method is only applied to model system with few degrees of freedom.

B. Ehrenfest mean-field method

A simple way of modeling non-adiabatic phenomena is the trajectory-based Ehrenfest mean-field (EMF) approach [32-34]. In this approach, the nuclear motion is governed by weighted averages of potential energy surfaces for the related BO states, and the time dependent weights in the average are determined by integrating electronic Schrödinger equation. For a system of $N$ atoms and $n$ electrons, the nuclei follow trajectory ${\bf R}(t)$=$({\bf R}_1, {\bf R}_2, \ldots{\bf R}_N)$, while the electronic wave function $\Phi({\bf r}, t)$ with $\bf r$=$({\bf r}_1, {\bf r}_2, \ldots{\bf r}_n)$ satisfies Schrödinger equation

 $i\hbar \frac{\partial }{\partial t}\Phi \left({\bf r}, t \right) = {\hat H_{\textrm{el}}}\left( \bf{r, R} \right)\Phi \left({\bf r}, t \right)$ (24)

Here, the partial electronic Hamiltonian is defined for fixed nuclear coordinate $\bf R$ as

 ${\hat H_{\textrm{el}}}\left( \bf{r, R} \right) = - \frac{{{\hbar ^2}}}{{2{m_\textrm{e}}}}\sum\limits_{i = 1}^n \nabla _i^2 + \hat V\left( \bf{r, R} \right)$ (25)

Again, the electronic wave function can be written as a linear combination of adiabatic eigen functions $\phi_j(\bf{r, R})$, i.e.,

 $\Phi \left({\bf r}, t \right) = \mathop \sum \limits_j {C_j}\left( t \right){\phi _j}\left( \bf{r, R} \right)$ (26)

Here, ${\phi _j}\left( \bf{r, R} \right)$ is the solution of the time-independent electronic Schrödinger equation

 ${\hat H_{\textrm{el}}}\left( \bf{r, R} \right){\phi _j}\left( \bf{r, R} \right) = {E_j}\left(\bf R \right){\phi _j}\left(\bf {r, R} \right)$ (27)

for fixed nuclear position $\bf R$ at time $t$. (Please note that the adiabatic state is denoted by index $j$ or $k$ rather than $\alpha$ or $\beta$ in this and the following subsections in order to be consistent with the literature.) Substituting the ansatz into the time dependent electronic Schrödinger equation followed by multiplication of $\phi^*_k(\bf{r, R})$ and integration over $\bf r$ leads to a set of coupled differential equation of $C_k$

 $i\hbar {\dot C_k} = \sum\limits_j {C_j}\left( t \right)\left( {{V_{kj}} - i\hbar \dot{\bf R} \cdot {{\bf d}_{kj}}} \right)$ (28)

with ${{\bf d}_{kj}}$=$\langle{\phi _k}\left( \bf{r, R} \right)\left|\displaystyle\frac{\partial }{{\partial\bf R}}\right|{\phi _j}\left( \bf{r, R} \right)\rangle$. The square modulus of the expansion coefficient $|C_k |^2$ is the probability of finding the system in adiabatic state $k$ at time $t$. It can be seen from Eq.(28) that the transitions between electronic states are promoted by the off-diagonal elements of Hamiltonian $V_{kj}$ and the non-adiabatic coupling ${\bf{\dot R}} \cdot {{\bf{d}}_{kj}}$.

To incorporate electronic-nuclear feedback, the atoms are assumed to be evolved on an effective potential representing an average over the adiabatic states weighted by their state population $|C_k |^2$, i.e.

 ${E_{{\rm{eff}}}} = \langle\Phi |{\hat H_{\textrm{el}}}|\Phi\rangle = \sum\limits_k |{C_k}{|^2}{E_k}$ (29)

This is the well-known Ehrenfest mean-field approach. Then the nuclear forces can be obtained from the gradient of the mean energy

 ${\bf F} = - \left\langle\Phi \left|\frac{\partial }{{\partial\bf R}}{\hat H_{\textrm{el}}}\right|\Phi\right\rangle\\~~~~~ = - \sum\limits_k |{C_k}{|^2}\frac{{\partial {E_k}}}{{\partial\bf R}} + \sum\limits_{k, j} C_k^*{C_j}\left( {{E_j} - {E_k}} \right){{\bf d}_{jk}}$ (30)

The first term is simply the average forces over the adiabatic states, and the second term involves non-adiabatic coupling between the adiabatic states.

Since the nuclear motion is represented by one point in phase space for all the electronic states in the EMF method, the mean-field trajectory can be dramatically different from the behavior of the quantum wave-packets if the potential energy surfaces have very different slopes for the related BO states. In addition, the effective potentials for forward and backward pathways are different in the Ehrenfest scheme, resulting in different populations, which violates microscopic [5].

C. Trajectory surface hopping

Unlike the EMF method, the nuclei move always on the potential energy surface of a single adiabatic state $\phi_k$ and switching from one state to another occurs suddenly in infinitesimal time in the trajectory surface-hopping (TSH) [39, 40]. The hopping is determined stochastically in a way to maintain consistency between the fraction of trajectories associated with a particular electronic state and the population on that state. A "fewest switch" surface hopping (FSSH) algorithm was firstly proposed to determine the transition among different electronic states [40]. Instead of calculating the state-averaged path in the EMF method, the TSH technique involves an ensemble of $N$ trajectories. The ensemble-averaged number of trajectory on adiabatic state $k$ at any time is equal to its occupation number, i.e.

 ${N_k}\left( t \right) = {\rho _{kk}}\left( t \right)N$ (31)

where the density matrix is defined by

 ${\rho _{kj}} = C_k^*\left( t \right){C_j}\left( t \right)$ (32)

At a short time $t$+$\delta t$ later, the new occupation number is

 ${N_k}\left( {t'} \right) = {\rho _{kk}}\left( {t'} \right)N$ (33)

By assuming $N_k(t')$$<$$N_k(t)$, then the minimum number of transition from $N_k(t)$ to $N_k(t')$ is $\delta N_k$=$N_k(t')$$-$$N_k(t)$ hops from state $\phi_k$ to any other state and zero hops from any other state to state $\phi_k$. The transition probability during this time interval $\delta t$ is

 ${P_k}\left( {t, \delta t} \right) = \frac{{\delta {N_k}}}{{{N_k}}} = \frac{{{\rho _{kk}}\left( t \right) - {\rho _{kk}}\left( {t'} \right)}}{{{\rho _{kk}}\left( t \right)}} = - \frac{{{{\dot \rho }_{kk}}\delta t}}{{{\rho _{kk}}\left( t \right)}}$ (34)

From Eqs. (28) and (32), we can obtain

 ${\dot \rho _{kj}} = - \sum\limits_l \bigg[{\frac{i}{\hbar }\left( {{{\rm{V}}_{kl}}{\rho _{lj}}-{\rho _{kl}}{V_{lj}}} \right) + } \\~~~~~~ {}{{\bf{\dot R}} \cdot \left( {{{\bf d}_{kl}}{\rho _{lj}}-{\rho _{kl}}{{\bf d}_{lj}}} \right)}\bigg]$ (35)

Then

 ${\dot \rho _{kk}}= \sum\limits_{l \ne k} {b_{kl}}\\~~~~~ = \sum\limits_{l \ne k} \left[{\frac{2}{\hbar }{\rm{Im}}\left( {\rho _{kl}^*{V_{kl}}} \right)-2{\rm{Re}}\left( {\rho _{kl}^*{\bf \dot R} \cdot {{\bf d}_{kl}}} \right)} \right]$ (36)

Since the probability $P_k(t, \delta t)$ must be the sum over the probability of each electronic state, $P_{kl}$, for a transition from $\phi_k$ to a specific state of $\phi_l$, thus

 ${P_k}\left( {t, \delta t} \right) = \sum\limits_l {P_{kl}}\left( {t, \delta t} \right)$ (37)

It follows that

 ${P_{kl}}\left( {t, \delta t} \right) = - \frac{{\displaystyle\frac{2}{\hbar }{\rm{Im}}\left( {\rho _{kl}^*{V_{kl}}} \right) - 2{\rm{Re}}\left( {\rho _{kl}^* {\bf\dot R} \cdot {{\bf d}_{kl}}} \right)}}{{{\rho _{kk}}}}$ (38)

A transition from $\phi_k$ to $\phi_l$ occurs if

 $P_k^{\left( {m - 1} \right)} < \varsigma < P_k^{\left( m \right)}$ (39)

where $\varsigma$ (0$\leq$$\varsigma$$\leq$1) is a uniform random number

and

 $P_k^{\left( m \right)} = \sum\limits_l^m {P_{kl}}s$ (40)

is the sum of transition probability for the first $m$ states. If switch occurs, the trajectory will now begin to evolve on the potential energy surface $V_{ll}(\bf R)$. At the same time, a velocity adjustment must be made at a position of transition $\bf R$ in order to conserve the total energy,

 ${{\bf P}_{\left( k \right)}} = {{\bf P}_{\left( l \right)}} + \Delta {\bf P} \cdot {{\bf d}_{kl}}$ (41)
 $\frac{{{\bf P}_{\left( k \right)}^2}}{{2M}} + {V_{kk}}\left( {\bf R} \right) = \frac{{{\bf P}_{\left( l \right)}^2}}{{2M}} + {V_{ll}}\left( {\bf R} \right)$ (42)

If a hop is not energetically allowed, the hopping attempt is ignored.

We took the trans-cis photoisomerization of the ethylene-bridged arobenzene as an example to perform EMF and TSH non-adiabatoc dynamics simulations [82]. The EMF equation of motion is solved with the nuclear motion numerically integrated with the velocity Verlet algorithm. Meanwhile, the laser pulse is characterized by a vector potential, which is coupled to the electronic Hamiltonian through the time dependent Peierls substitution [82]. The combined electronic structure calculations and non-adiabatic dynamics simulations reveal that the initial relaxation from the S$_1$ state of the trans isomer is the downhill motion to the region of conical intersection (CI) with a lifetime of $\sim$30 fs, which is consistent with the ultrafast S$_1$ deactivation reported experimentally. However, the S$_1$ lifetime ($\sim$1.6 ps) for the cis isomer was overestimated, as compared with the first time constant of $\sim$70 fs inferred experimentally.

Subsequently, a systematic study on the photoisomerization dynamics of the ethylene-bridged arobenzene in the gas phase and the CH$_3$OH solution was performed by the combined quantum mechanics with surface-hopping (QM-SH) and molecular mechanics (QM-SH/MM) [84]. As shown in FIG. 2, the QM moiety of the ethylene-bridged arobenzene was treated with ab initio based surface-hopping method, while the MM moiety of the CH$_3$OH molecules was simulated with the OPLS-AA force field. The interaction between the QM and MM moieties was described by electrostatic-embedding scheme with the total energy of the entire system given by

 ${E_{{\rm{tot}}}} = {E_{{\rm{QM}}}} + {E_{{\rm{MM}}}} + E_{{\rm{QM}}/{\rm{MM}}}^{{\rm{el}} + {\rm{vdW}}}$ (43)
 FIG. 2 The QM/MM scheme.

where $E_{\rm{QM/MM}}^{\rm{el+vdW}}$ contains the contribution from electrostatic and van der Waals (vdW) interactions between the QM and MM moieties. The effective Hamiltonian of the QM moiety taking into account the MM environment can be written as

 ${\hat H_{{\rm{eff}}}} = {\hat H_0} + {\hat H'}\\~~~~ {\hat H'} ={\hat {V'}_{Nn}} + {\hat {V'}_{Ne}}\\~~~~~~~~ = \sum\limits_{N \in {\rm{QM}}}\sum\limits_{n \in {\rm{MM}}} \frac{{{q_n}{Z_n}}}{{{{\bf R}_n} - {{\bf R}_N}}} -\\~~~~~~~~~~~~ {}\sum\limits_{i \in {\rm{QM}}\left( {{\rm{electrons}}} \right)} \sum\limits_{n \in {\rm{MM}}} \frac{{{q_n}}}{{{{\bf R}_n} - {{\bf r}_i}}}$ (44)

where $\hat H_0$ is the standard QM Hamiltonian, $\hat H'$ is the electrostatic interaction of nuclear charges $Z_N$ and electrons $i$ in QM with the MM point charges $q_n$. ${\bf R}_n$ is the position vector of point charge in the MM moiety, while ${\bf R}_N$ and ${\bf r}_i$ are position vectors of nuclei and electrons in the QM region, respectively. The effective Hamiltonian was used to calculate energies, energy gradients, and non-adiabatic coupling vectors for the hybrid QM/MM surface hopping implementation.

On the basis of the QM-SH/MM dynamics simulation, the S$_1$ lifetime was fitted to be 31.2 fs for the trans isomer in the gas phase, which is nearly unchanged by the solvent effect. However, the S$_1$ dynamics of the cis isomer is considerably influenced by the solvent effect with the S$_1$ lifetime increased from 80.6 fs in the gas phase to 127.4 fs in the CH$_3$OH solution. After the S$_1$$\rightarrowS_0 transition, isomerization to the cis isomer in the gas phase has a time constant of 510 fs, which is decreased to 420 fs by the solvent effect. The time scales for isomerization to the trans isomer in the gas phase and the CH_3OH solution are similar to each other with a value of \sim300 fs. The QM-SH/MM simulated results indicate that the solvent effect would improve the back-and-forth photo-switch efficiency between the cis and trans isomers and the trans\rightarrowcis or cis\rightarrowtrans photoisomerization has a high quantum yield (0.4) [82]. Therefore, the ethylene-bridged arobenzene may have good photo-responsive properties as a photoactive unit in real systems. D. Decoherence correction for the MQC non-adiabatic dynamics simulation In a full quantum description, decoherence occurs when an initial wave packet is divided into a couple of wave packets moving on separate adiabatic states. Because the spawned wave packets experience different potentials, the nuclear dynamics for different quantum pathway diverge in both position and phase. The random interaction leads to destructive interference between the nuclear wave functions that are associated with these pathways. However, the nuclear motion is represented by one point in phase space for all the electronic states in the EMF method, the mean-field trajectory can be dramatically different from the behavior of the quantum wave-packets if the potential energy surfaces have very different slopes for the related BO states. Like the EMF method, the nuclear motion is also treated classically and the electronic wave functions of different BO states retain coherence over the entire course in each trajectory in TSH method. The coherent effect is overestimated in the TSH and EMF algorithms. To solve this problem, various decoherence corrections have been developed, which are briefly described as follows. A straightforward approach considering decoherence correction in the MQC non-adiabatic dynamics methods is to introduce a damping term which collapses the coherent wave function onto a single adiabatic state. Tully proposed a simple decoherence correction scheme [40] by damping the off diagonal elements of the density matrix of Eq.(35) with  \dot \rho _{kj}^D = - \gamma \left( {1 - {\delta _{kj}}} \right){\rho _{kj}} (45) Here parameter \gamma is the decoherence rate. In Tully's original approach, the damping term is arbitrarily parameter-dependent, which needs to be properly determined. A proper treatment of quantum coherence in non-adiabatic dynamics requires explicit attention to the decoherence time scale. For time longer than the decoherence interval, the memory of the phase relationship among the quantum states is lost. Therefore, the decoherence time \tau_k can be defined by [62]  \frac{1}{{{\tau _k}}} = \sum\limits_{j \ne k}^N {\left| {{C_j}\left( t \right)} \right|^2}{\gamma _{kj}} (46) Here, the sum is over all the adiabatic state except the state under consideration. The decoherence rates for each pair of state, \gamma_{kj}, have been approximately determined by various methods. One of them is based on the semi-classical approximation for the nuclear degree of freedom [41, 42]. Supposing the nuclear wave function being the superposition of the frozen Gaussian, Schwartz and coworkers calculated the nuclear overlap integral, and defined the quantum decoherence rate as  \frac{1}{{{\tau _{kj}}}} = {\gamma _{kj}} = \sum\limits_n \frac{{\left| {{{\bf F}_k}\left( t \right) - {{\bf F}_j}\left( t \right)} \right|}}{{2\hbar \sqrt {{\alpha _n}} }} (47) Here, {\bf F}_k(t)$$-$${\bf F}_j(t) is the Hellmann-Feynman force difference between two classical trajectories, and \alpha_n is the Gaussian width parameter for the nth nuclear coordinates. In Ref.[42],  {\alpha _n} = \frac{{6{M_n}{k_{\rm{B}}}T}}{{{\hbar ^2}}} (48) is chosen with M_n being the effective mass of the nth particle. This expression is only valid at high temperature for displaced harmonic oscillator. By using the self-consistent decay-of-mixing semi-classical trajectory method [35, 36], Zhu and Truhlar gave the decoherence time as  {\tau _{kj}} = \left| {\frac{\hbar }{{\Delta {E_{kj}}}}\left( {1 + \frac{C}{{{E_{{\rm{kin}}}}}}} \right)} \right| Here, \Delta E_{kj} is the potential energy difference between two states in the diabatic representation, E_{\rm{kin}} is the kinetic energy of the nuclei, and C is an empirical parameter. Alternatively, by calculating the electronic density matrix based on the Gaussian wave packet anzatz, Shenvi, Subotinik and Yang [57, 61] derived a decoherence rate as  {\gamma _{kj}} = {\rm{Re}}\left( {{\alpha _{kj}}\Delta {x_{kj}}\Delta {{\dot x}_{kj}}} \right) (50) Here, \alpha_{kj}=\alpha_k\alpha_J^*/(\alpha_k+\alpha_j^*) with \alpha_k being the width parameter of the Gaussian wave, \Delta x_{kj} is the spatial separation between Gaussians, and \Delta {\dot x_{kj}} is the rate of separation between Gaussian centers. Starting from quantum-classical Liouville equation, Subotnik and coworkers derived the decoherence rate for a 2-level system [56, 65, 66] as  \begin{array}{l} \gamma _{12}^{\left( 1 \right)} \approx \displaystyle\frac{1}{2} \sum\limits_n \left( {{\bf F}_{11}^n - {\bf F}_{22}^n} \right)\displaystyle\frac{1}{{\rho _{12}^{\left( 1 \right)}}}\frac{{\partial \rho _{12}^{\left( 1 \right)}}}{{\partial {{\bf P}_n}}}\\ \gamma _{12}^{\left( 2 \right)} \approx - \displaystyle\frac{1}{2}\sum\limits_n \left( {{\bf F}_{11}^n - {\bf F}_{22}^n} \right)\frac{1}{{\rho _{12}^{\left( 2 \right)}}}\displaystyle\frac{{\partial \rho _{12}^{\left( 2 \right)}}}{{\partial {{\bf P}_n}}} \end{array} (51) Here, {\bf F}_{ij}^n is the general Hellmann-Feynman forces for the nth nuclear mode, \rho_{ij}^{(\lambda)} is the density matrix elements for the \lambdath active electronic state, and {\bf P}_n is the momentum of the nth nuclear mode. Different algorithms have been implemented to deal with decoherence in the TSH simulations [62]. For example, the decoherence time for state \phi_k is computed either from t_0 or the time of the previous decoherence event in the decoherence-induce surface hopping approach [62]. If the current time step t_m is within the coherence interval for each basis state, a random number with unit of \tau_k is generated. If the number is larger than the time interval between the current time and the time of the previous decoherence event, the state vector is reduced. A uniformly distributed random number is generated and compared to the quantum probability |C_k|^2 of the adiabatic state \phi_k at time t_m. If the random number is less than or equal to the quantum probability, the nuclear trajectory hops to the decohering state \phi_k, and the wave function is reset to \phi_k. Otherwise, \phi_k is projected out from the current state vector, and the wave function is renormalized. Meanwhile, the nuclear velocities are reset as in standard TSH algorithm. To take into account the quantum decoherence effect without calculating the decoherence time, Granucci and coworkers set up an overlap decoherence correction (ODP) method [54]. In addition to the representative phase space point {\bf Q}_j^n, {\bf P}_j^n traveling on the current PES j for a given trajectory n, ancillary points {\bf Q}_{k, I}^n, {\bf P}_{k, I}^n are attributed to every other electronic state k. By assigning frozen Gaussian wave packets to the representative point, the electronic population is modified by evaluating the overlap between wave packets k, I, and j, which depends on their phase space distance. If the overlap drops below a given threshold, then the wave packet is suppressed and the coefficient of current state is replaced. When a surface hopping occurs, a new leading Gaussian wave packet is created, and all the other wave packets on the same potential energy surface are suppressed. This method depends on two parameters, the width of the wave packet and the threshold overlap, and is meaningful only for short time and near degeneration region. A recent application of the ODP method is the photoisomerization of arylazopyrazole photoswitches [83]. The simulations predicted a new stereospecific unidirectional excited-state relaxation mechanism. Stochastic evolution of the wave function in surface hopping creates additional computational challenge. To solve this limitation, mean-field approaches accounting for decoherence effects have been developed. Schwartz and his coworkers proposed a mean-field dynamics approach with stochastic decoherence (MF-SD) [37]. As in traditional Ehrenfest dynamics, the initial density matrix elements are propagated by Eq.(28), and the mean-field Hellmann-Feynman forces Eq.(30) are acting on each nucleus. Then the decoherence probability for each adiabatic state is calculated by  {P_k} = \frac{{{\rho _{kk}}}}{{{\tau _k}}}{\rm d}t (52) Like surface hopping dynamics, collapse of the system is determined by comparing the probability to a random number. If collapse is allowed, the mean-field HF force is reset to adiabatic force for state \phi_k, the electronic density matrix is reduced to a pure state, and the classical velocity is rescaled. This algorithm implies that MF-SD is not a fully deterministic MF method. Another MF method with decoherence is the so-called coherent switching method with decay-of-mixing (CSDM) introduced by Zhu et al. [35, 36]. In their method, density matrix is propagated with two terms. One comes from fully coherent contribution as standard Ehrenfest term, and the other is the decay-of-mixing term which is similar to Eq.(45), where the decoherence time (or rate) is determined by Eq.(49). The equations of motion of the nuclear positions are also governed by two terms, the coherent mean force of Eq.(30) and the decoherent force that concerns the angular momentum of the mapped electronic variable  {{\bf F}^D} = - \frac{{M{{\dot V}^D}}}{{{\bf P} \cdot {\bf\hat s}}}{\bf\hat {s}} (53) Here, M is the mass, \dot V^D is the time derivative of the decoherent effective potential energy, \bf P is the momentum, and \bf\hat s is the unit vector in that direction energy which is deposited or consumed. Though the decaying of the off-diagonal elements are deterministic in the CSDM method, stochastic surface hops, similarly to FSSH, are then undertaken. Therefore, the CSDM is a SH-MF combined method, which still requires averaging over a large number of stochastic realizations. Recently, Prezhdo's group developed a fully deterministic MF method with decoherence, the so-called coherence penalty functional (CPF) method [38]. To account for decoherence effects, they introduced an argument Hamiltonian  {\tilde H_{\textrm{el}}} = {\hat H_{\textrm{el}}} + \sum\limits_{i \ne j} {\lambda _{kj}}{\left| {C_k^*{C_j}} \right|^2} (54) The additional term penalizes the coherence development. The penalty between pair of state i and j is determined by a constant \lambda_{kj}, which is proportional to the decoherence rate \gamma_{kj}=1/\tau_{kj}  {\lambda _{kj}} = f\frac{1}{{{\tau _{kj}}}} = f\frac{{\left| {{{\bf F}_k}\left( t \right) - {{\bf F}_j}\left( t \right)} \right|}}{{2\hbar \sqrt {{\alpha _0}} }} (55) Empirically, f=100/P_0 and \alpha_0=1 Bohr are chosen. Here P_0 is the initial nuclear momentum. The coherence penalty term contributes a partial correction \lambda_{kj}\rho_{kj} to the mean-field dynamics. Then standard Ehrenfest techniques are performed. No spatial redistribution or state collapse of the wave function is included in CPF method, therefore, it is simple in implementation and computationally efficient. E. Quantum trajectory mean-field Method The aforementioned decoherence corrections typically enter the equations through coupling terms involving parameters that characterize the nuclear relaxation process. We recently developed a quantum trajectory mean-field method (QTMF) to account for the dephasing effects in a parameter-free fashion. Since the classical nuclear motion can be regarded as continuum quantum measurement to the electronic subsystem, the feedback of nuclear motion on the motion of electron should be involved at any time. According to this idea, we expanded the electronic wave function in the following format,  \psi \left( {{\bf{R, r}}, t} \right) = \sum\limits_j {C_j}\left( {{{\bf R}_t}, \;t} \right)|{\phi _j}\left( {{{\bf R}_t}, {\bf r}} \right)\rangle (56) Here, the coefficient {C_j}\left( {{{\bf R}_t}, \;t} \right) not only explicitly depends on time t, but also explicitly on the nuclear potion {\bf R}_t={\bf R}(t) at time t. This condition implies that C_j is different if the nuclear position changes. We substituted this wave function into the electronic Schrödinger Eq.(24) to yield the differential equation of {C_k}\left( {{{\bf R}_t}, \;t} \right),  i\hbar \frac{{\partial {C_k}\left( {{{\bf R}_t}, \;t} \right)}}{{\partial t}} = {C_k}\left( {{{\bf R}_t}, \;t} \right){E_k}\left( {{{\bf R}_t}} \right) -\\~~~~ i\hbar\sum\limits_jC_j({\bf R}_t, t){\bf\dot R}\cdot {\bf d}_{kj}-i\hbar {\nabla _R}{C_k}\left( {{{\bf R}_t}\;t} \right) \cdot {\bf\dot R}\quad\quad (57) Compared to Eq.(28), this time evolution equation has an additional term contributed by the evolution of classical degree of freedom. This term was also obtained in the Liouville equation [26]. As a matter of fact, if we consider the expansion coefficient C_j({\bf R}_t, {\bf P}_t, t) to be depending explicitly both on the nuclear position and momentum, we can recover the Liouville equation. Nevertheless, in the framework of mixed quantum-classical trajectory method, the nuclear momentum is uniquely determined by its position, so we assume that C_j(R_t, t) is only (R_t, t)-dependent. When the interactions between the adiabatic states are assumed weak, Eq.(57) can be solved iteratively as a perturbation expansion. To the zeroth order perturbation, i.e., when the coupling between two states are ignored, the formal solution of C_k({\bf R}_t, t) is  \begin{array}{l} {C_k}\left( {{{\bf{R}}_t},t} \right) = {C_k}\left( {{{\bf{R}}_{t - \delta t}},t - \delta t} \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left[ { - \frac{i}{\hbar }\int_{t - \delta t}^t {\rm{d}} t'{E_k}\left( {{{\bf{R}}_{t'}}} \right)} \right] \end{array} (58) Substituting this solution back into Eq.(58), we obtain  \begin{array}{l} {\nabla _R}{C_k}\left( {{{\bf{R}}_t},t} \right) = {C_k}\left( {{{\bf{R}}_{t - \delta t}},t - \delta t} \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\exp \left[ { - \frac{i}{\hbar }\int_{t - \delta t}^t {\rm{d}} t'{E_k}\left( {{{\bf{R}}_{t'}}} \right)} \right] \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( { - \frac{i}{\hbar }} \right){\nabla _R}\left( {\int_{t - \delta t}^t {\rm{d}} t'{E_k}\left( {{{\bf{R}}_t}} \right)} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = {C_k}\left( {{{\bf{R}}_t},t} \right) \cdot \left( { - \frac{i}{\hbar }} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\nabla _R}\left( {\int_{t - \delta t}^t {\rm{d}} t'{E_k}\left( {{{\bf{R}}_t}} \right)} \right) \end{array} (59) By using the mean value theorem of integrals, we yield  \int_{t - \delta t}^t {\rm{d}}t'{E_k}\left({\bf R}_t \right) = {E_k}\left( {\bf\bar R} \right)\delta t (60) where \bf\bar R is the particular nuclear coordinate in the time interval \delta t. In the limit of small time interval, E_j(\bf\bar R) approximately equals to E_j({\bf R}_t). Finally, we derive the evolution equation of C_k({\bf R}_t, t) as  \begin{array}{l} i\hbar \frac{{\partial {C_k}\left( {{{\bf{R}}_t},t} \right)}}{{\partial t}} = {C_k}\left( {{{\bf{R}}_t},t} \right){E_k}\left( {{{\bf{R}}_t}} \right) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i\hbar \sum\limits_j {{C_j}} \left( {{{\bf{R}}_t},t} \right){\bf{\dot R}} \cdot {{\bf{d}}_{kj}} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\Delta {\bf{R}} \cdot {\nabla _R}{E_k}\left( {{{\bf{R}}_t}} \right){C_k}\left( {{{\bf{R}}_t},t} \right) \end{array} (61) In the format of density matrix with \rho _{kj}\left({\bf R}_t, t\right)= C_k^*\left( {\bf R}_t, t \right){C_j}\left({\bf R}_t, t\right), it reads  \dot \rho _{kj} = - \frac{i}{\hbar }\left( E_k\left( {\bf R}_t \right) - E_j\left( {\bf R}_t \right) \right)\rho _{kj} -\\~~~~~~~~~~~~~~ {} \sum\limits_l \dot{\bf R} \cdot \left( {\bf d}_{kl}\rho _{lj} - \rho _{kl}{\bf d}_{lj} \right)+\\~~~~~~~~~~~~~~ {}\frac{i}{\hbar }\Delta{\bf R} \cdot \left(\nabla _RE_k\left( {\bf R}_t \right) - \nabla _RE_j\left( {\bf R}_t \right) \right)\rho _{kj} (62) In the last term, \Delta\bf R is the nuclear displacement in one time step \delta t, \nabla_RE_k({\bf R}_t)$$-$$\nabla_RE_j({\bf R}_t) is the Hellmann-Feynman force difference between two adiabatic states, and their product is the work of nuclear movement acting on the electrons. This work well describes the influence of the classical trajectories back onto the electronic state evolution, we call it back-action work. We will see that this term not only contributes to the phase of the electronic wave function, but also induces decoherence. In the mean-field trajectory method, the forces on the nuclei (see Eq.(30)) involve both the average forces over the adiabatic states and the non-adiabatic coupled forces. This implies that nuclear motion continuously extracts information from the electronic state, including the energy and phase information. Due to the coupling between the adiabatic states, the mean force has both real and imaginary parts, so does the nuclear displacement \Delta\bf R. These facts imply that the real part of classical back-action work in Eq.(61) delays the electronic phase, while its imaginary part induces collapse of the electronic superposition state. Moreover, this dephasing term is free of any parameters and can be easily implemented in simulations. We applied our method to calculate the probability of transmission of three well-known Tully's model systems, i.e., the single-avoided crossing model, the dual-avoided crossing model, and the extended coupling model. The initial samples were prepared in the lower states by uniformly select 2000 position points. At each position, proper momentum was set so that the 2000 samples have the same energy. The four-order Range-Kuttle algorithm was used for the propagation of the electronic wave function (Eq.(62)) with integration step of 0.1 fs. For the integration of the nuclear motion (Eq.(30)), velocity-Verlet method was employed. The results are given in FIG. 3. For the single and dual-avoided crossing models, our method gave results in good agreement with FSSH and exact methods. It indicates that the back-action work of nucleus on the electron well describes the dephasing effects. The advantage of our method is that the dephasing term is derived directly from the electronic Schrödinger equation, and is free of any parameters. Unfortunately, the present method fails for the extended coupling model; no reflection was observed. We infer this failure as a consequence of zeroth order perturbation approximation. For such strongly coupled system, higher order expansion terms in the formal solution of C_k(R_t, t) should be involved. A proper solution is the Lindblad-type master equation.  FIG. 3 The probability of transmission of (a) single-avoided crossing model, (b) dual-avoided crossing model, and (c) extended coupling model. From an idea of quantum measurement picture, the classical nuclei can be viewed as classical observers to the electronic subsystem. Based on this picture, a Lindblad-type master equation was derived for non-adiabatic MD simulation [69],  \dot \rho \left( t \right) = - \frac{i}{\hbar }\left[{\tilde H, \rho \left( t \right)} \right] + \sum\limits_{j \ne k} {\gamma _{F, jk}}{\cal D}\left[{{M_{jk}}\left( R \right)} \right]\rho \left( t \right)+\\~~~~~~~~~~~~~~ {}\sum\limits_{j \ne k} \sqrt {{\gamma _{F, jk}}} {\cal H}\left[{{M_{jk}}\left( R \right)} \right]\rho \left( t \right){\xi _{jk}}\left( t \right) (63) Here \tilde H is the effective Hamiltonian with element  \tilde H_{kj} = {E_k}\delta _{kj} - i\hbar {\bf\dot R} \cdot {\bf d}_{kj} (64) The Lindblad super-operator {\cal D} has form  {\cal D}\left[M_{jk}(R)\right]\rho(t)= M_{jk}\rho (t)M_{jk}^\dagger \ -\frac{1}{2}\left\{ M_{jk}^\dagger \ M_{jk}, \rho(t) \right\}\\ (65) with  {M_{jk}}\left( R \right) = |{\phi _j}\left({\bf R} \right)\rangle\langle{\phi _j}\left({\bf R} \right)| - |{\phi _k}\left({\bf R} \right)\rangle\langle{\phi _k}\left({\bf R} \right)| (66) The {\cal H} operator is of form  {\cal H}\left[{{M_{jk}}\left( R \right)} \right]\rho \left( t \right) = {M_{jk}}\rho \left( t \right) + \rho \left( t \right){M_{jk}}^\dagger - \\~~~~~~~~~~~~~~~~~~~~~~~~~~ {}{\rm{Tr}}\left[{\left({M_{jk}}^\dagger + {M_{jk}}\right)\rho \left( t \right)} \right] (67) These two operators represent the information gain that the nuclei motion back-acting on the electrons. They lead to decoherence between states \phi_j(\bf R) and \phi_k(\bf R). Since the nuclei motion is much slower than its electronic counterpart, a coarse grained fluctuating force associating with adiabatic potential energy surface is supposed to evolve the Newtonian equation. By assuming the average force satisfying a Gaussian distribution, the information gain rate coinciding with decoherence rate is obtained, which is expressed as  \gamma _{F, kj} = \frac{{{{\left( {{{\bf F}_k}\left( t \right) - {{\bf F}_j}\left( t \right)} \right)}^2}}}{{{{\left( {|{{\bf F}_k}\left( t \right)| + \left| {{{\bf F}_j}\left( t \right)} \right|} \right)}^2}}}\frac{1}{{{\tau _c}}} (68) Here the characteristic time \tau_c scales atomic motion, which is typically of picosecond. With respect to the femtosecond timescale of the electronic part, in practice we adopt 1/\tau_\textrm{c}$$\approx$10$^{-3}$ in (reduced) units of the electronic energies. The density matrix dynamics also include an additional term that accounts for the quantum-jump related stochastic noise, $\xi_{kj}(t)$. This term has a role of collapsing the electronic state from a quantum superposition onto a single basis state. The information gain rate given by Eq.(68) is of configuration dependence, but it does not need any microscopic information of the nuclear (quantum) wave packets. It allows thus for a convenient implementation even for simulation on complex molecular systems. This protocol provides a natural interface between the separate quantum and classical treatment. Even for the demanding extended coupling model, it gives satisfactory reflection probabilities [69].

The QTMF approach is numerically implemented on the fly with ab initio calculation at the level of the complete active space self-consistent field, which is employed to explore photodissociation dynamics of diazirinone [84]. For comparison, the photodissociation process has been simulated with the TSH method and ab initio multiple spawning (AIMS) method. Upon irradiation of N$_2$CO at $\sim$355 nm, the dissociation was predicted to proceed in a stepwise way. The first C$-$N bond fission occurs along the S$_1$ pathway, while the second C$-$N bond is broken in the S$_0$ state as a result of the S$_1$$\rightarrowS_0 transition in the vicinity of the first conical intersection (CI-1). The time constants for cleavage of the first and second C-N bonds were predicted to be respectively 73.0 and 137.0 fs by the QTMF simulations, which are close to the AIMS simulated values of 67.3 and 112.3 fs respectively. In comparison, the time constants for fission of the C-N bonds were significantly overestimated by the TSH simulation. From the AIMS simulated wave functions, the transition ratio was deduced to be 5.2% for the S_1$$\rightarrow$S$_0$ transition in the region of the second conical intersection (CI-2), which is in good agreement with the corresponding value of 4.6% obtained by the QTMF simulation. However, the transition ratio is estimated to be 8.3% by the TSH simulation. As discussed before, the coherent effect is mainly responsible for the overestimation of the time constants for the C$-$N fissions and the transition ratio by the TSH simulation. Since the decoherence correction is treated quantum mechanically in the AIMS simulation, a good agreement in the QTMF and AIMS simulated results provides a clear evidence that direct ab initio QTMF approach is a reliable tool for simulating nonadiabatic dynamics processes [84].

Schrödinger equation is replaced with a master equation in the QTMF scheme [84], where the quantum decoherence is automatically fulfilled by the stochastic process. Therefore, the QTMF method is more suitable for simulating non-adiabatic processes that involve more than two BO states, as compared with the TSH approach. Experimentally, it has been established that at least the lowest three singlet states are involved in photodynamics of acetylacetone (AcAc) at $\sim$265 nm. Therefore, we took AcAc as a representative system to carry out direct ab initio QTMF non-adiabatic dynamics simulation [85]. Upon irradiation at $\sim$265 nm, the S$_2$ state is initially populated in the Franck-Condon region. Three types of non-adiabatic processes were determined by the combined ab initio calculation and QTMF simulation, which are summarized in FIG. 4.

 FIG. 4 Three types of non-adiabatic processes upon irradiation of the AcAc to the S$_2$ state: (a) S$_2$($^1\pi\pi^*$)$\rightarrow$S$_2$-1$\rightarrow$ S$_2$/S$_1$-1$\rightarrow$S$_1$$\rightarrowS_1/S_0$$\rightarrow$S$_0$, (b) S$_2$($^1\pi\pi^*$)$\rightarrow$S$_2$-2$\rightarrow$S$_2$/S$_1$-2$\rightarrow$S$_1$$\rightarrowS_1/S_0$$\rightarrow$S$_0$, and (c) S$_2$($^1\pi\pi^*$)$\rightarrow$S$_2$/S$_1$/S$_0$$\rightarrowS_0. There is a strong coherence between the S_2 and S_1 states in the initial in-plane proton transfer FIG. 4(a) and the out-of-plane rotational isomerization FIG. 4 (b) and (c). The subsequent decoherence takes place mainly in the S_2/S_1-1 intersection region, which is the predominant route to the S_1 state. Meanwhile, the S_1 state is populated through the S_2$$\rightarrow$S$_1$ transition in the conical intersection S$_2$/S$_1$-2 region. Although the chelated enolic structure is the most stable in the S$_1$ state, the chelated structure produced in the S$_1$ state is gradually changed into the non-chelated isomers. The ratio of the chelated/non-chelated isomers is about 2:3 in the S$_1$ state. A similar ratio is found in the S$_0$ state. Formation of the non-chelated isomer in the S$_0$ state is considerably overestimated by the present QTMF simulation, as compared with the experimental findings ($\sim$2:1). The first time component was experimentally inferred to be $\sim$50.0 fs in the gas phase or in dioxane, acetonitrile, and $n$-hexane, which is well reproduced by the QTMF simulation. The related process is exclusively assigned as the in-plane proton transfer from the S$_2$ Franck-Condon structure to the local minimum (E-S$_2$-1). However, a considerable difference in the second time constant of the S$_2$ state was reported in the previous experimental studies. The underlying process for the second time constant was determined to be the rotational isomerization to the global minimum (E-S$_2$-2) by the QTMF simulation. This rotamerization was characterized to proceed in the S$_2$ and S$_1$ superposition state, which is probably the reason for the difference reported experimentally for the second time constant of the S$_2$ state. The S$_1$ lifetime was estimated to be 2.11 ps by the present ab initio QTMF simulation, which is in excellent agreement with the experimentally inferred value of 2.12, 2.13, and 2.25 ps for acetylacetone in $n$-hexane, acetonitrile, and dioxane respectively [85].

F. Factorization of electron-nuclear wave function

The recently developed coupled-trajectory mixed quantum-classical (CT-MQC) approach by Gross's group [86] and non-Hermitian surface hopping (nH-SH) method suggested by Gao and Thiel [68] are also able to treat electronic decoherence in non-adiabatic processes. These two methods based on Born-Huang expansion of the total wave function, in which the nuclear wave function is written in polar form with amplitude and phase part. The modulus of the nuclear wave function is constructed in terms of frozen Gaussian-shape wave packets, while the phase is connected to the classical momentum. Decoherence corrections containing a term referring to quantum momentum, $p_{n, I}^{\textrm{qm}}$=$-i\hbar\nabla_n\chi_I/\chi_I$ (here $\chi_I$ is the nuclear wave function at $I$th state and its derivative is for the $n$th nuclear coordinate), are introduced in both the electronic density matrix evolution equation and nuclear Newtonian equation. In CT-MQC and nH-SH, the quantum momentum is replace by its classical approximation $p_{n, I}^{\textrm{qm}}$=$-\alpha_n[{\bf R}(t)$$-$${\bf R}_I (t)]$ with $\alpha_n$ being the width of the Gaussian wave packet. The CT-MQC method employs mean-filed algorithm to solve the equations. The difference of CT-MQC approach from other mean field approaches is that the calculation of quantum momentum needs to collect information about the positions of all trajectories at a given time. The nH-SH method adopts fewest switch algorithm. Like QTMF method, the factorization of the electron-nuclear wave function in these methods catches the main features of the time-dependent potential and can correctly describe the nuclear dynamic, but the evolvement of a combination of Gaussian-shaped wave packets makes their numerical implementation more complicated than QTMF.

G. Multiple spawning method

To describe the wave function splitting, Martinez and coworkers developed an approximate ab initio multiple spawning (AIMS) method [70-75], which is founded from a quantum mechanical standpoint, yet retains a strong resemblance to classical mechanism. The time dependent wave function is represented as

 $\Psi =\sum\limits_I {\phi _I}\left( \bf{r, R} \right){\chi _I}\left( {\bf R}, t \right)$ (69)

The wave function of electronic states $\phi_I(\bf{r, R})$ are taken to be orthonormal over the electronic coordinates

 $\int{\rm{d}}r\phi _J^*\left(\bf{r, R} \right){\phi _I}\left( \bf{r, R} \right) = {\delta _{IJ}}$ (70)

The time dependent nuclear wave function for each electronic state $I$, $\chi_I({\bf R}, t)$, is expanded in MS approach as a linear combination of localized travelling Gaussian wave functions

 ${\chi _I}\left({\bf R}, t\right) = \sum\limits_j C_j^I\left( t \right)\chi _j^I\left( {{\bf R};{\bf\bar R}_j^I\left( t \right), {\bf\bar P}_j^I\left( t \right), \bar \gamma _j^I\left( t \right), \bar \alpha _j^I} \right)\\$ (71)

Here, the indices $j$ and $I$ label the $j$th nuclear basis function on the $I$th electronic state. Each of the Gaussian wave function is constructed as a product of one-dimensional Gaussian basis functions

 $\chi _j^I\left( {{\bf R};{\bf\bar R}_j^I\left( t \right), {\bf\bar P}_j^I\left( t \right), \bar \gamma _j^I\left( t \right), \bar \alpha _j^I} \right) =\\~~~~~~~~~~~~ {\rm{e}}^{i\bar \gamma _j^I( t )} \sum\limits_{\lambda = 1}^{3N} \tilde \chi _j^I\left(R_\lambda;\bar R_{\lambda j}^I( t), \bar P_{\lambda j}^I( t ), \bar \alpha _{\lambda j}^I \right)$ (72)

with

 $\tilde \chi _j^I\left( {{R_\lambda };\bar R_{\lambda j}^I\left( t \right), \bar P_{\lambda j}^I\left( t \right), \bar \alpha _{\lambda j}^I} \right)={\left( {\frac{{2\bar \alpha _{\lambda j}^I}}{\pi }} \right)^{1/4}}\cdot\\~~~~ \exp \left[{-\bar \alpha _{\lambda j}^I{{\left( {{R_\lambda }-\bar R_{\lambda j}^I\left( t \right)} \right)}^2} + i\bar P_{\lambda j}^I\left( t \right)\left( {{R_\lambda }-\bar R_{\lambda j}^I\left( t \right)} \right)} \right]\\$ (73)

Here, the index $\lambda$ represents the Cartesian coordinates of the nuclei ($x$, $y$, and $z$). Each Gaussian basis function is parametrized with time dependent position $\bar R_{\lambda j}^I(t)$, momentum $\bar P_{\lambda j}^I(t)$, nuclear phase $\bar\gamma_{\lambda j}^I(t)$, and a time independent width $\bar\alpha_{\lambda j}^I$. The equations of motion of the parameters are given by

 $\begin{array}{l} \displaystyle\frac{{\partial \bar R_{\lambda j}^I\left( t \right)}}{{\partial t}} = \frac{{\bar P_{\lambda j}^I\left( t \right)}}{{{M_\lambda }}}\\ \displaystyle\frac{{\partial \bar P_{\lambda j}^I\left( t \right)}}{{\partial t}} = - {\left[{\frac{{\partial {V_{II}}\left({\bf R} \right)}}{{\partial R}}} \right]_{\bar R_{\lambda j}^I}}\\ \displaystyle\frac{{\partial \bar \gamma _j^I\left( t \right)}}{{\partial t}} = - {V_{II}}\left( {{\bf\bar R}_j^I\left( t \right)} \right) + \sum\limits_{\lambda = 1}^{3N} \frac{{{{\left( {\bar P_{\lambda j}^I\left( t \right)} \right)}^2} - 2\bar \alpha _{\lambda j}^I}}{{2{M_\lambda }}} \end{array}$ (74)

This is equivalent to an interaction picture with a semi-classical reference. The time evolution of the complex coefficient is governed by the nuclear Schrödinger equation

 $\frac{{{\rm{d}}C_j^I\left( t \right)}}{{{\rm{d}}t}} = - i\sum\limits_{kl} {\left( {{\bf S}_{II}^{ - 1}} \right)_{kl}} \times\\~~~~~~~~~~~~ \left[{{{\left( {{{\bf H}_{II}}-i{{\bf\dot S}_{II}}} \right)}_{k, l}}C_l^J + \sum\limits_{J \ne I} {{({{\bf H}_{IJ}})}_{kl}}C_l^J} \right]$ (75)

Here, ${{\bf S}_{II}}$=$\langle\chi _k^I{\rm{|}}\chi _l^I\rangle$ and ${\bf\dot S}_{II}$=$\left\langle {\chi _k^I|(\partial /\partial t\chi _l^I)} \right\rangle$ are the time-dependent nuclear overlap matrix and its time derivative, and $({\bf H}_{IJ})_{kl}$=$\left\langle {\chi _k^I|{\bf\hat H}|\chi _l^J} \right\rangle$ is the Hamiltonian matrix describing the interstate coupling between basis function on electronic state $I$ and $J$. The electronic amplitudes determine the population $|C_I|^2$, and the degree of coherence of different electronic states, i.e., the elements of reduced density matrix ${\rho _{KJ}}$=$C_K^*\left( t \right){C_J}\left( t \right)$.

In principle, the basis set consists of an infinite nuclear basis functions. However, such implementation is impractical. To provide a procedure which is quasi-classical, an ensemble of initial conditions is generated such that the coordinate and momentum space distributions of the ensemble can mimic the quantum mechanical initial state. In practice, "virtual" basis functions that are close to the trajectory of the initial state in phase space during the non-adiabatic event are chosen. For initial conditions, usually only one electronic state is populated. The initial position and momentum that define each of the Gaussian basis functions are sampled from the Wigner distribution. On the other hand, a basis set of fixed size may not well describe the quantum mechanical phenomena. A practical implementation of the spawning method is to control the growth of the basis set by allowing it to expand only when the dynamics signals impending failure of classical mechanics, such as non-adiabatic or tunneling effects. This can be achieved by monitoring the magnitude of the non-adiabatic coupling for each nuclear basis function. When the nuclear basis function enters strong non-adiabatic coupling region, the solution of nuclear Schrödinger equation stops and new basis functions are spawned (or generated), with zero initial population on the other electronic state. The position and momentum of the newly spawned basis functions are determined by using Frank-Condon considerations, i.e., they have maximal overlap with their parent basis function at that time point. The spawning method sews classical flesh onto quantum mechanical bones and gives concrete meaning to multistate dynamics. AIMS has been successfully applied to study the photo-induced isomerization of various molecules [70-73, 87].

In the spirit of AIMS, Shalashilin's group develops the so-called ab initio multiple cloning (AIMC) algorithm [75] for quantum non-adiabatic molecular dynamics based on multiple configurational Ehrenfest (MCE) method [74]. Instead of using classical dynamics to guide the basis coherent states as in AIMS, MCE uses Ehrenfest mean-field trajectories with the Hamiltonian averaged over two or more electronic PESs to guide the frozen Gaussian basis functions. To overcome the difficulty of standard Ehrenfest method that only includes a single configuration and is not suitable for accurate non-adiabatic dynamics, MCE method expands the exact molecular wave function as a sum of multiple wave packets (or configurations). Each Ehrenfest configuration includes electronic and nuclear part. The nuclear motion is described by a single frozen Gaussian wave packet as in AIMS method. Cloning procedure is introduced when any of the trajectory basis functions (TBFs) become strongly mixed, i.e., when a TBF passes a conical intersection or region of strong non-adiabatic coupling. The cloning does not change the wave function, it only re-expresses the mixed states Ehrenfest configuration as a superposition of distinct TBFs, each is dominated by a single electronic state. The AIMC method can be viewed as a commentary approach to AIMS method.

Ⅲ. SUMMARY

We have present an overview of the recent progress in the theoretical treatment of electronically non-adiabatic processes by using mixed quantum-classical methods. The accuracy and efficiency of these methodologies can be summarized as follows. Although full and accurate quantum dynamics theories are powerful algorithms to propagate nuclear wave packets, these techniques require that the Hamiltonian can be expanded as a sum of product of one-particle operator, or the single-particle function basis can be transformed to discrete variable representation basis. For systems larger than six atoms, it is not computationally affordable to perform full quantum dynamical calculation. A proper approximation is to treat the electrons or protons as quantum part and the other nuclear degrees of freedom by classical mechanics. Following from a partial Wigner transform over the classical degrees of freedom of the Liouville equation for the full quantum system, mixed quantum-classical Liouville equation correctly accounting for the coupled evolution of the quantum subsystem and classical subsystem is obtained. However, the numerical solution of this equation needs to sample a huge number of trajectories in the phase space, which limits its application only to a few degrees of freedom.

The Ehrenfest mean-field and trajectory surface hopping algorithms are the two most widely used MQC approaches for the treatment of non-adiabatic dynamics. In these methods, the motions of classical particle produce time variations of the quantum particles' Hamiltonian and induce transitions among the quantum states. These changes of quantum states, in turn, alter the forces that govern the classical motion. With mean-field method, trajectories evolve on a averaged potential energy surface, while the surface-hopping trajectories evolve on a single potential energy surface, with instantaneous hops between surfaces. Because the classical trajectories do not account for the wave packet separation, the phase in EMF and TSH methods are over coherent. A number of procedures including decoherence have been introduced. However, the decoherence effects typically enter equations through coupling terms involving parameters that characterized the bath relaxation process. The recently developed quantum trajectory mean-field method by us partially solve this problem. It is appropriate for weakly coupled system.

The mean-field and surface hopping methods try to sew quantum mechanical flesh onto classical bones. Alternatively, factorization of the electron-nuclear wave function, such as coupled-trajectory mixed quantum-classical, non-Hermitian surface hopping, and multiple spawning methods attempt to sew classical flesh onto quantum mechanical bones. They expand the nuclear wave function for an electronic state in terms of Gaussian wave packets, while the propagation of the wave packets is guided by classical dynamics. These methods provide promising approaches to both large systems and non-adiabatic dynamics.

Beside dealing with the non-adiabaticity as reviewed in this article, theoretical treatment of non-adiabatic dynamics also needs to solve the multistate electronic structure problem. Appropriate ab initio multi-reference quantum chemical method, such as generalized valence bond (GVB), multi-configuration self-consistent field (MCSCF), complete active space SCF (CASSCF), full configuration interaction (FCI) theories are usually employed for the evaluation of both the gradient and non-adiabatic coupling for the excited states. More efficient time-dependent density functional theory (TDDFT) and multi-configuration pair-density functional theory (MC-PDFT) have also been developed for strongly correlated system. The electronic structure theories go beyond this article, we refer to Ref.[88] for a more comprehensive review.

In conclusion, the non-adiabatic dynamics theories and practical implementation of numerical methods have received a significant boost in the past decades, but still leave a lot of room for future improvements.

Ⅳ. ACKNOWLEDGEMENTS

This work was supported by the National Key R & D Program of China (No.2017YFB0203405) and the National Natural Science Foundation of China (No.21421003).

Reference
 [1] M. Born, and K. Huang, Dynamical Theory of Crystal Lattices. Oxford, UK: Oxford University Press (1954). [2] D. M. Hirst, Potential Energy Surfaces: Molecular Structure and Reaction Dynamics. London: Taylor & Francis (1985). [3] M. Baer, Beyond Born-Oppenheimer: Electronic Nona-diabatic Coupling Terms and Conical Intersections. Hoboken, New Jersey: John Wiley & Son (2006). [4] A. W. Jasper, C. Y. Zhu, S. Nangiaa, and D. G. Truhlar, Faraday Discuss. 127 , 1 (2004). DOI:10.1039/b405601a [5] J. C. Tully, J. Chem. Phys. 137 , 22A301 (2012). [6] G. A. Worth, and L. S. Cederbaum, Annu. Rev. Phys. Chem. 55 , 127 (2004). DOI:10.1146/annurev.physchem.55.091602.094335 [7] T. Yonehara, K. Hanasaki, and K. Takatsuka, Chem. Rev. 112 , 499 (2012). DOI:10.1021/cr200096s [8] S. Deb, and P. M. Weber, Annu. Rev. Phys. Chem. 62 , 19 (2011). DOI:10.1146/annurev.physchem.012809.103350 [9] W. Domcke, and D. R. Yarkony, Annu. Rev. Phys. Chem. 63 , 325 (2012). DOI:10.1146/annurev-physchem-032210-103522 [10] J. H. Lehman, and M. I. Lester, Annu. Rev. Phys. Chem. 65 , 537 (2014). DOI:10.1146/annurev-physchem-040513-103628 [11] W. J. Schreier, P. Gilch, and W. Zinth, Annu. Rev. Phys. Chem. 66 , 497 (2015). DOI:10.1146/annurev-physchem-040214-121821 [12] H. S. Goan, and G. J. Milburn, Phys. Rev. B 64 , 235307 (2001). DOI:10.1103/PhysRevB.64.235307 [13] L. D. Landau, Phys. Z. Sowjetunion 2 , 46 (1932). [14] C. Zener, Proc. R. Soc. A 137 , 696 (1932). DOI:10.1098/rspa.1932.0165 [15] M. Topaler, and N. Makri, J. Chem. Phys. 100 , 4430 (1996). DOI:10.1021/jp951673k [16] R. B. Gerber, V. Buch, and M. A. Ratner, J. Chem. Phys. 77 , 3022 (1982). [17] R. Kosloff, J. Chem. Phys. 92 , 2087 (1988). DOI:10.1021/j100319a003 [18] P. Jungwirth, and R. B. Gerber, J. Chem. Phys. 102 , 8855 (1995). DOI:10.1063/1.468939 [19] P. Jungwirth, and R. B. Gerber, J. Chem. Phys. 102 , 6046 (1995). DOI:10.1063/1.469339 [20] P. Jungwirth, E. Fredj, and R. B. Gerber, J. Chem. Phys. 104 , 9332 (1996). DOI:10.1063/1.471678 [21] M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324 , 1 (2000). DOI:10.1016/S0370-1573(99)00047-2 [22] H. Meyer, and G. A. Worth, Theor. Chem. Acc. 109 , 251 (2003). DOI:10.1007/s00214-003-0439-1 [23] H. B. Wang, and M. Thoss, J. Chem. Phys. 119 , 1289 (2003). DOI:10.1063/1.1580111 [24] D. Hochstuhl, C. M. Hinz, and M. Bonitz, Eur. Phys. J. Spec. Top. 223 , 177 (2014). DOI:10.1140/epjst/e2014-02092-3 [25] H. B. Wang, J. Phys. Chem. A 119 , 7951 (2015). DOI:10.1021/acs.jpca.5b03256 [26] R. Kapral, and G. Ciccotti, J. Chem. Phys. 110 , 8919 (1999). DOI:10.1063/1.478811 [27] S. Nielsen, R. Kapral, and G. Ciccotti, J. Chem. Phys. 112 , 6543 (2000). DOI:10.1063/1.481225 [28] M. Santer, U. Manthe, and G. Stock, J. Chem. Phys. 114 , 2001 (2001). DOI:10.1063/1.1336576 [29] K. Ando, Chem. Phys. Lett. 360 , 240 (2002). DOI:10.1016/S0009-2614(02)00848-5 [30] I. Horenko, C. Salzmann, B. Schmidt, and C. Schütte, J. Chem. Phys. 117 , 11075 (2002). DOI:10.1063/1.1522712 [31] K. Ando, and M. Santer, J. Chem. Phys. 118 , 10399 (2003). DOI:10.1063/1.1574015 [32] J. W. Negele, Rev. Mod. Phys. 54 , 913 (1982). DOI:10.1103/RevModPhys.54.913 [33] G. D. Billing, Chem. Phys. Lett. 100 , 535 (1983). DOI:10.1016/0009-2614(83)87423-5 [34] N. L. Doltsinis, and D. Marx, J. Theor. Comput. Chem. 1 , 319 (2002). DOI:10.1142/S0219633602000257 [35] C. Y. Zhu, S. Nangia, A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. 121 , 7658 (2004). DOI:10.1063/1.1793991 [36] C. Y. Zhu, A. W. Jasper, and D. G. Truhlar, J. Chem. Theory Comput. 1 , 527 (2005). DOI:10.1021/ct050021p [37] M. J. Bedard-Hearn, R. E. Larsen, and B. J. Schwartz, J. Chem. Phys. 123 , 234106 (2005). DOI:10.1063/1.2131056 [38] A. V. Akimov, R. Long, and O. V. Prezhdo, J. Chem. Phys. 140 , 194107 (2014). DOI:10.1063/1.4875702 [39] J. C. Tully, and R. K. Preston, J. Chem. Phys. 55 , 562 (1971). DOI:10.1063/1.1675788 [40] J. C. Tully, J. Chem. Phys. 93 , 1061 (1990). DOI:10.1063/1.459170 [41] E. Neria, and A. Nitzan, J. Chem. Phys. 99 , 1109 (1993). DOI:10.1063/1.465409 [42] B. J. Schwartz, E. R. Bittner, O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 104 , 5942 (1996). DOI:10.1063/1.471326 [43] O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 107 , 5863 (1997). DOI:10.1063/1.474312 [44] O. V. Prezhdo, and P. J. Rossky, J. Chem. Phys. 107 , 825 (1997). DOI:10.1063/1.474382 [45] Y. L. Volobuev, M. D. Hack, M. S. Topaler, and D. G. Truhlar, J. Chem. Phys. 112 , 9716 (2000). DOI:10.1063/1.481609 [46] M. D. Hack, and D. G. Truhlar, J. Chem. Phys. 114 , 2894 (2001). DOI:10.1063/1.1342224 [47] C. Y. Zhu, K. Nobusada, and H. Nakamura, J. Chem. Phys. 115 , 3031 (2001). DOI:10.1063/1.1386811 [48] C. Y. Zhu, H. Kamisaka, and H. Nakamura, J. Chem. Phys. 115 , 11036 (2001). DOI:10.1063/1.1421070 [49] C. Y. Zhu, H. Kamisaka, and H. Nakamura, J. Chem. Phys. 116 , 3234 (2002). DOI:10.1063/1.1446032 [50] I. R. Craig, and D. E. Manolopoulos, J. Chem. Phys. 121 , 3368 (2004). [51] C. Y. Zhu, A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. 120 , 5543 (2004). DOI:10.1063/1.1648306 [52] A. W. Jasper, and D. G. Truhlar, J. Chem. Phys. 123 , 064103 (2005). DOI:10.1063/1.1995695 [53] J. P. Rank, and R. Kapral, J. Chem. Phys. 132 , 074106 (2010). DOI:10.1063/1.3310811 [54] G. Granucci, M. Persico, and A. Zoccante, J. Chem. Phys. 133 , 134111 (2010). DOI:10.1063/1.3489004 [55] J. E. Subotnik, J. Phys. Chem. A 115 , 12083 (2011). DOI:10.1021/jp206557h [56] J. E. Subotnik, and N. Shenvi, J. Chem. Phys. 134 , 024105 (2011). DOI:10.1063/1.3506779 [57] N. Shenvi, J. E. Subotnik, and W. T. Yang, J. Chem. Phys. 134 , 144102 (2011). DOI:10.1063/1.3575588 [58] P. Shushkov, R. Li, and J. C. Tully, J. Chem. Phys. 137 , 22A549 (2012). DOI:10.1063/1.4766449 [59] S. Fernandez-Alberti, A. E. Roitberg, T. Nelson, and S. Tretiak, J. Chem. Phys. 137 , 014512 (2012). DOI:10.1063/1.4732536 [60] B. R. Landry, and J. E. Subotnik, J. Chem. Phys. 137 , 22A513 (2012). DOI:10.1063/1.4733675 [61] N. Shenvi, and W. T. Yang, J. Chem. Phys. 137 , 22A528 (2012). [62] H. M. Jaeger, S. Fischer, and O. V. Prezhdo, J. Chem. Phys. 137 , 22A545 (2012). DOI:10.1063/1.4757100 [63] A. Kelly, and T. E. Markland, J. Chem. Phys. 139 , 014104 (2013). DOI:10.1063/1.4812355 [64] V. N. Gorshkov, S. Tretiak, and D. Mozyrsky, Nat. Commun. 4 , 2144 (2013). [65] J. E. Subotnik, W. J. Ouyang, and B. R. Landry, J. Chem. Phys. 139 , 214107 (2013). DOI:10.1063/1.4829856 [66] J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. J. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67 , 387 (2016). DOI:10.1146/annurev-physchem-040215-112245 [67] C. Y. Zhu, Sci. Rep. 6 , 24198 (2016). DOI:10.1038/srep24198 [68] X. Gao, and W. Thiel, Phys. Rev. E 95 , 013308 (2017). DOI:10.1103/PhysRevE.95.013308 [69] W. Feng, L. T. Xu, X. Q. Li, W. H. Fang, and Y. J. Yan, AIP Adv. 4 , 077131 (2014). DOI:10.1063/1.4891821 [70] T. J. Martinez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100 , 7884 (1996). DOI:10.1021/jp953105a [71] M. Ben-Nun, and T. J. Martnez, J. Chem. Phys. (108), 7244 (1998). [72] M. Ben-Nun, J. Quenneville, and T. J. Martnez, J. Phys. Chem. A 104 , 5161 (2000). [73] B. G. Levine, and T. J. Martnez, Annu. Rev. Phys. Chem. 58 , 613 (2007). DOI:10.1146/annurev.physchem.57.032905.104612 [74] K. Saita, and D. V. Shalashilin, J. Chem. Phys. 137 , 22A506 (2012). DOI:10.1063/1.4734313 [75] D. V. Makhov, W. J. Glover, T. J. Martinez, and D. V. Shalashilin, J. Chem. Phys. 141 , 054110 (2014). DOI:10.1063/1.4891530 [76] H. D. Meyera, and W. H. Miller, J. Chem. Phys. 70 , 3214 (1979). DOI:10.1063/1.437910 [77] G. Stock, and M. Thoss, Phys. Rev. Lett. 78 , 578 (1997). DOI:10.1103/PhysRevLett.78.578 [78] W. H. Miller, J. Phys. Chem. A 105 , 2942 (2001). DOI:10.1021/jp003712k [79] W. H. Miller, J. Phys. Chem. A 113 , 1405 (2009). DOI:10.1021/jp809907p [80] W. H. Miller, and S. J. Cotton, J. Chem. Phys. 145 , 081102 (2016). DOI:10.1063/1.4961551 [81] J. Liu, J. Chem. Phys. 145 , 204105 (2016). DOI:10.1063/1.4967815 [82] L. H. Liu, S. Yuan, W. H. Fang, and Y. Zhang, J. Phys. Chem. A 115 , 10027 (2011). [83] Y. T. Wang, X. Y. Liu, G. L. Cui, W. H. Fang, and W. Thiel, Angew. Chem. Int. Ed. Engl. 128 , 14215 (2016). DOI:10.1002/ange.v128.45 [84] B. B. Xie, L. H. Liu, G. L. Cui, W. H. Fang, J. Cao, W. Feng, and X. Q. Li, J. Chem. Phys. 143 , 194107 (2015). DOI:10.1063/1.4935800 [85] B. Xie, G. Cui, and W. Fang, J. Chem. Theory Comput. 13 , 2717 (2017). DOI:10.1021/acs.jctc.7b00153 [86] S. K. Min, F. Agostini, and E. K. Gross, Phys. Rev. Lett. 115 , 73001 (2015). DOI:10.1103/PhysRevLett.115.073001 [87] L. Liu, S. Xia, and W. Fang, J. Phys. Chem. A 118 , 8977 (2014). DOI:10.1021/jp5019923 [88] L. Gagliardi, D. G. Truhlar, G. L. Manni, R. K. Carlson, C. E. Hoyer, and J. L. Bao, Acc. Chem. Res. 50 , 66 (2017). DOI:10.1021/acs.accounts.6b00471