The article information
 Lianghui Gao, Binbin Xie, Weihai Fang
 高靓辉, 谢斌斌, 方维海
 Theories and Applications of Mixed QuantumClassical Nonadiabatic Dynamics
 混合量子经典非绝热动力学理论与应用
 Chinese Journal of Chemical Physics, 2018, 31(1): 1226
 化学物理学报, 2018, 31(1): 1226
 http://dx.doi.org/10.1063/16740068/31/cjcp1712234

Article history
 Received on: December 6, 2017
 Accepted on: January 6, 2018
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 socalled BornOppenheimer (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 fulldynamics simulations can be carried onthefly 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 [36]. 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 nonadiabatic processes are also essential parts of collisions of electronically excited species, chemiluminescent reactions, recombination reactions, heterolytic dissociations, and electron transfer processes [611]. Nonadiabatic dynamics also play intimate roles in quantum information processing with atomic, molecular, or solidstate qubits [12]. Proper incorporation of nonadiabatic 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 nonadiabatic 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 twostates 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
$ 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) 
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 singlesurface dynamics occurs in the region of configuration spaces where the two PESs nearly cross. Under the assumption that near the conical crossing point (
$ {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,
The LZformula provides a good approximation for nonadiabaticity 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 nonadiabatic 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 nonadiabatic 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 quantumclassical (MQC) approaches where the electrons are treated as quantum particle while the nuclei are treated as classical particles, and semiclassical method that use a semiclassical description for all the degree of freedom of electrons and nuclei. Timedependent Hartree and multiconfiguration timedependent Hartree are fully quantum mechanical methods that propagate wave packets for nuclei [1625]. 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 nonadiabatic 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 quantummechanical description, while the heavy classical particles are represented by a trajectory description. Mixed quantumclassical Liouville (QCL) equation in the form of density matrix is a promising method [2631]. A variety of quantumclassical trajectory approaches have also been proposed, including Ehrenfest meanfield (EMF) method [3238], trajectory surface hopping (TSH) method [3968], quantumtrajectory meanfield (QTMF) method [69], and multiple spawning (MS) method [7075]. In this review, we summarize the most significant development of the mixed quantumclassical methods and their applications. Because the coupling of quantum and classical DOF in mixed quantumclassical approaches is not completely consistent, an alternative strategy is to use a semiclassical description for all the degrees of freedom [7681]. As proposed by Meyer and Miller, the electronic state is characterized by a pair of classic harmonic oscillator actionangle variable, which can be transformed into Cartesian coordinates and momenta for the actual numerical trajectory calculation. We suggest the readers to refer to Refs.[7681] for a detailed discussion and review of semiclassical method.
Ⅱ. NONADIABATIC DYNAMICS THEORIES AND APPLICATIONS A. Mixed quantumclassical Liouville equationFor 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 quantumclassical system is composed of an
$ \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 Qcoordinates is essentially performed on the density matrix
$ {\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
$ \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
$ \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 quantumclassical Liouville equation.
To solve the evolution equation of
$ {\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}\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
$ \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 nonadiabatic 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 HellmannFeynman 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
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}_{\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 firstorder 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
The QCL approach provides a correct quantumclassical 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 meanfield methodA simple way of modeling nonadiabatic phenomena is the trajectorybased Ehrenfest meanfield (EMF) approach [3234]. 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
$ 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
$ {\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 \left({\bf r}, t \right) = \mathop \sum \limits_j {C_j}\left( t \right){\phi _j}\left( \bf{r, R} \right) $  (26) 
Here,
$ {\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
$ 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
To incorporate electronicnuclear feedback, the atoms are assumed to be evolved on an effective potential representing an average over the adiabatic states weighted by their state population
$ {E_{{\rm{eff}}}} = \langle\Phi {\hat H_{\textrm{el}}}\Phi\rangle = \sum\limits_k {C_k}{^2}{E_k} $  (29) 
This is the wellknown Ehrenfest meanfield 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 nonadiabatic 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 meanfield trajectory can be dramatically different from the behavior of the quantum wavepackets 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 hoppingUnlike the EMF method, the nuclei move always on the potential energy surface of a single adiabatic state
$ {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
$ {N_k}\left( {t'} \right) = {\rho _{kk}}\left( {t'} \right)N $  (33) 
By assuming
$ {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}\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
$ P_k^{\left( {m  1} \right)} < \varsigma < P_k^{\left( m \right)} $  (39) 
where
and
$ P_k^{\left( m \right)} = \sum\limits_l^m {P_{kl}}s $  (40) 
is the sum of transition probability for the first
$ {{\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 transcis photoisomerization of the ethylenebridged arobenzene as an example to perform EMF and TSH nonadiabatoc 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 nonadiabatic dynamics simulations reveal that the initial relaxation from the S
Subsequently, a systematic study on the photoisomerization dynamics of the ethylenebridged arobenzene in the gas phase and the CH
$ {E_{{\rm{tot}}}} = {E_{{\rm{QM}}}} + {E_{{\rm{MM}}}} + E_{{\rm{QM}}/{\rm{MM}}}^{{\rm{el}} + {\rm{vdW}}} $  (43) 
where
$ {\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
On the basis of the QMSH/MM dynamics simulation, the S
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 meanfield trajectory can be dramatically different from the behavior of the quantum wavepackets 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 nonadiabatic 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
A proper treatment of quantum coherence in nonadiabatic 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
$ \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,
$ \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,
$ {\alpha _n} = \frac{{6{M_n}{k_{\rm{B}}}T}}{{{\hbar ^2}}} $  (48) 
is chosen with
$ {\tau _{kj}} = \left {\frac{\hbar }{{\Delta {E_{kj}}}}\left( {1 + \frac{C}{{{E_{{\rm{kin}}}}}}} \right)} \right $ 
Here,
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,
Starting from quantumclassical Liouville equation, Subotnik and coworkers derived the decoherence rate for a 2level 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,
Different algorithms have been implemented to deal with decoherence in the TSH simulations [62]. For example, the decoherence time for state
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
Stochastic evolution of the wave function in surface hopping creates additional computational challenge. To solve this limitation, meanfield approaches accounting for decoherence effects have been developed. Schwartz and his coworkers proposed a meanfield dynamics approach with stochastic decoherence (MFSD) [37]. As in traditional Ehrenfest dynamics, the initial density matrix elements are propagated by Eq.(28), and the meanfield HellmannFeynman 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 meanfield HF force is reset to adiabatic force for state
Another MF method with decoherence is the socalled coherent switching method with decayofmixing (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 decayofmixing 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,
Recently, Prezhdo's group developed a fully deterministic MF method with decoherence, the socalled 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
$ {\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,
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 meanfield method (QTMF) to account for the dephasing effects in a parameterfree 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
$ 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
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
$ \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
$ \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
$ \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,
In the meanfield trajectory method, the forces on the nuclei (see Eq.(30)) involve both the average forces over the adiabatic states and the nonadiabatic 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
We applied our method to calculate the probability of transmission of three wellknown Tully's model systems, i.e., the singleavoided crossing model, the dualavoided 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 fourorder RangeKuttle 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)), velocityVerlet method was employed. The results are given in FIG. 3. For the single and dualavoided crossing models, our method gave results in good agreement with FSSH and exact methods. It indicates that the backaction 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
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 Lindbladtype master equation was derived for nonadiabatic 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_{kj} = {E_k}\delta _{kj}  i\hbar {\bf\dot R} \cdot {\bf d}_{kj} $  (64) 
The Lindblad superoperator
$ {\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}\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 backacting on the electrons. They lead to decoherence between states
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
The QTMF approach is numerically implemented on the fly with ab initio calculation at the level of the complete active space selfconsistent 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
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 nonadiabatic 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
There is a strong coherence between the S
The recently developed coupledtrajectory mixed quantumclassical (CTMQC) approach by Gross's group [86] and nonHermitian surface hopping (nHSH) method suggested by Gao and Thiel [68] are also able to treat electronic decoherence in nonadiabatic processes. These two methods based on BornHuang 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 Gaussianshape wave packets, while the phase is connected to the classical momentum. Decoherence corrections containing a term referring to quantum momentum,
To describe the wave function splitting, Martinez and coworkers developed an approximate ab initio multiple spawning (AIMS) method [7075], 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
$ \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
$ {\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
$ \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
$ \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 semiclassical 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,
In principle, the basis set consists of an infinite nuclear basis functions. However, such implementation is impractical. To provide a procedure which is quasiclassical, 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 nonadiabatic 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 nonadiabatic or tunneling effects. This can be achieved by monitoring the magnitude of the nonadiabatic coupling for each nuclear basis function. When the nuclear basis function enters strong nonadiabatic 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 FrankCondon 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 photoinduced isomerization of various molecules [7073, 87].
In the spirit of AIMS, Shalashilin's group develops the socalled ab initio multiple cloning (AIMC) algorithm [75] for quantum nonadiabatic 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 meanfield 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 nonadiabatic 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 nonadiabatic coupling. The cloning does not change the wave function, it only reexpresses 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.
Ⅲ. SUMMARYWe have present an overview of the recent progress in the theoretical treatment of electronically nonadiabatic processes by using mixed quantumclassical 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 oneparticle operator, or the singleparticle 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 quantumclassical 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 meanfield and trajectory surface hopping algorithms are the two most widely used MQC approaches for the treatment of nonadiabatic 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 meanfield method, trajectories evolve on a averaged potential energy surface, while the surfacehopping 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 meanfield method by us partially solve this problem. It is appropriate for weakly coupled system.
The meanfield and surface hopping methods try to sew quantum mechanical flesh onto classical bones. Alternatively, factorization of the electronnuclear wave function, such as coupledtrajectory mixed quantumclassical, nonHermitian 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 nonadiabatic dynamics.
Beside dealing with the nonadiabaticity as reviewed in this article, theoretical treatment of nonadiabatic dynamics also needs to solve the multistate electronic structure problem. Appropriate ab initio multireference quantum chemical method, such as generalized valence bond (GVB), multiconfiguration selfconsistent field (MCSCF), complete active space SCF (CASSCF), full configuration interaction (FCI) theories are usually employed for the evaluation of both the gradient and nonadiabatic coupling for the excited states. More efficient timedependent density functional theory (TDDFT) and multiconfiguration pairdensity functional theory (MCPDFT) 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 nonadiabatic 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.
Ⅳ. ACKNOWLEDGEMENTSThis work was supported by the National Key R & D Program of China (No.2017YFB0203405) and the National Natural Science Foundation of China (No.21421003).
[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 BornOppenheimer: Electronic Nonadiabatic 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/annurevphyschem032210103522 
[10]  J. H. Lehman, and M. I. Lester, Annu. Rev. Phys. Chem. 65 , 537 (2014). DOI:10.1146/annurevphyschem040513103628 
[11]  W. J. Schreier, P. Gilch, and W. Zinth, Annu. Rev. Phys. Chem. 66 , 497 (2015). DOI:10.1146/annurevphyschem040214121821 
[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/S03701573(99)000472 
[22]  H. Meyer, and G. A. Worth, Theor. Chem. Acc. 109 , 251 (2003). DOI:10.1007/s0021400304391 
[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/e2014020923 
[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/S00092614(02)008485 
[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/00092614(83)874235 
[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. BedardHearn, 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. FernandezAlberti, 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/annurevphyschem040215112245 
[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. BenNun, and R. D. Levine, J. Phys. Chem. 100 , 7884 (1996). DOI:10.1021/jp953105a 
[71]  M. BenNun, and T. J. Martnez, J. Chem. Phys. (108), 7244 (1998). 
[72]  M. BenNun, 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 