A multi-state trajectory method for non-adiabatic dynamics simulations

Coupled nuclear-electron motions are the central theme in chemical dynamics.^1–3 When the Born-Oppenheimer (BO) approximation breaks down, nonadiabatic transition occurs among different electronic states. The coupled nuclear-electron dynamics thus requires a multiple state description. A straightforward approach to include the contributions from different states is Ehrenfest dynamics⁴ in which nuclear motions are governed by the forces generated from the averaged potential energy surface (PES) of multiple electronic states. This mean field method may become invalid as the BO approximation recovers and the single state PES deviates appreciably from the averaged one. Surface hopping,^5–11 instead, incorporates electronic transitions among different states explicitly into the classical molecular dynamics. The transition between the coupled states is manipulated based on the presumed fewest switch algorithm to maintain the correct state population, and the transition probability is determined by the time-dependent electronic probability self-consistently. The overestimated electronic coherence in Ehrenfest dynamics is alleviated in surface hopping in which a swarm of classical trajectories travel along a single PES all the time except for fewest sudden switches among different electronic states. However the picture of the instantaneous state switch of a single classical nuclear trajectory in the early surface hopping treatment may not describe quantum coherence properly. A great deal of effort has been devoted to overcome the “decoherence” problem^12,13 in the surface hopping approach and the improvements could be substantial.¹¹ Recently Martens proposed a fully coherent stochastic surface hopping method in the quantum Liouville representation,¹⁴ which represents each density matrix element by an ensemble of trajectories to account for the quantum coherence.

In principle a full quantum description on both electronic and nuclear degrees of freedom (DOFs) in coupled nuclear-electron dynamics, such as the multiconfiguration time-dependent Hartree (MCTDH)¹⁵ and the multilayer MCTDH method,¹⁶ should be able to capture the quantum coherence exactly, and the requirement on the ad hoc hopping scheme may be therefore dismissed. However it is still a big challenge to treat nonadiabatic dynamics in realistic molecular systems in practice. To bypass the exponential scaling problem regarding to the grid basis size, the trajectory guided basis sets, such as coherent state, were suggested in multidimentional quantum dynamics simulations. For example, a number of promising trajectory based approaches have been developed, such as the matching-pursuit (MP)/split-operator-Fourier-transform (SOFT) method¹⁷ based on the adaptive coherent state expansion of the wavepacket, the multiple spawning approach^18,19 featured in the spawning of basis functions due to strong nonadiabatic couplings, and the approximate but less expensive multiconfigurational Ehrenfest (MCE) method²⁰ combining the coupled coherent states and Ehrenfest dynamics. At variance with the multiple spawning method, a path branching algorithm²¹ was suggested to describe quantum wavepacket bifurcation in nonadiabatic dynamics.

To account for quantum coherence accurately in a trajectory representation, it is important to treat the electronic and nuclear DOFs at the same dynamical footing. The Meyer-Miller (MM) model²² explains such a philosophical idea perfectly and treats the coupled nuclear-electron dynamics consistently by mapping the electronic DOFs into classical variables. Later on, Stock and Thoss^23,24 demonstrated that the MM Hamiltonian is the exact quantum representation for the coupled nuclear-electron system. The mapping scheme has been successfully implemented to a variety of nonadiabatic dynamical processes.^25–34 For example when the initial value representation (IVR) methodology^35,36 is applied to dynamics, the semiclassical description on the coupled nuclear-electron system has been demonstrated to be capable of reproducing quantum effects quite accurately. Depending on the level of approximations made to IVR, successful applications to interesting chemical or biochemical systems containing from tens to thousands of DOFs manifest the robustness of the methodology.^28–34 

As an extension of the equal footing idea, recently Cotton and Miller suggested a symmetrical windowing formalism^37,38 to treat the nuclear-electron dynamics quasiclassically. This symmetrical quasi-classical (SQC) approach retains the simplicity of the classical molecular dynamics (MD) simulation while incorporating electron transition inherently based on the MM model, therefore it is able to generate results in excellent agreement with exact quantum calculations on several benchmark model systems and appears very promising in the implementation to realistic complex molecular systems.^38–42 

To date, nuclear dynamics in nonadiabatic problems essentially evolves via one individual trajectory. In Ehrenfest dynamics, the single trajectory propagates on the averaged PES. Aside the instantaneous state switch, the individual nuclear trajectory in surface hopping is always on a particular state. The semiclassical or quasiclassical trajectory in the MM model is Ehrenfest-like, and the action-angle variables for the electronic DOFs determine which state the single trajectory is on. Of course, the final results are obtained based on an ensemble average of such individual trajectories.

In this work, we suggest a multi-state trajectory (MST) method to account for the fact that nuclei follow different dynamics on different states in nonadiabatic systems. For each state, there is one individual trajectory associated with it, which evolves only on the corresponding state, and the electron transition is incorporated by applying the MM model. The multi-state representation in the MST method allows each individual trajectory to propagate largely independently, therefore excellent numerical stability could be achieved especially for realistic molecular systems⁴³ in which the dynamics on different states may be so different that the coupled nuclear-electron dynamics could become quickly diverged or a very small time step has to be applied in other nonadiabatic simulation approaches such as SQC. For example, to model the electron coupled lithium ion diffusion in LiFePO₄, each state corresponds to the iron ion (Fe³⁺) locating in different sites, therefore the dynamics on different states differs substantially due to the strong Coulomb interactions, which causes the instability of the coupled nuclear-electron dynamics in the SQC simulation. Interesting comparison with the MCE approach²⁰ is noticed. The MCE approach introduces a multiconfiguration coherent state to describe nonadiabatic dynamics. However, each configuration is an Ehrenfest configuration, and the dynamics of the coherent state is also Ehrenfest. Our MST is also apparently in contrast to the hierarchical ones in the multiple spawning^18,19 or the path branching approaches.²¹ 

We derive the MST formalism in the MM framework and make approximations to obtain the quasiclassical version in Sec. II. The quasiclassical MST is applied in Sec. III to the spin-boson problem and two Tully’s models for nonadiabatic dynamics, and results and discussions are presented. Sec. IV concludes.

We start with the time-dependent Schrodinger equation for the coupled nuclear-electron system,

and the total wavefunction is written in terms of a product ansatz, i.e.,

with the φ_k(r) and χ_kl(R) to be the normalized electronic wavefunction and the lth normalized nuclear wavefunction for the state k, respectively. The index l denotes the labeling for the complete set of nuclear eigenfunctions, such as coherent states. Plugging Eq. (1b) into Eq. (1a) and multiplying it by $φ j ( r )$ results in the coupled differential equations for nuclear DOFs,

To obtain a classical-like dynamics, we write the nuclear wavefunction in terms of the action-angle variables, i.e.,

with A_kl and S_kl the amplitude and the phase, respectively. Substituting Eq. (2) into Eq. (1c), we obtain the following equation of motions (EOMs) isomorphic to the classical Hamilton-Jacobi equation

Introducing the transformation $P j l ≡∇ S j l R$⁠, we therefore arrive to the classical-like EOM, i.e.,

If the individual trajectory picture is taken, i.e., $χ j l R = χ k l R$⁠, Eq. (4) reduces to the well known Ehrenfest dynamics. And if $χ j l R ∣ χ k l R = δ j k$⁠, Eq. (4) represents for the BO dynamics. The MST approach assumes that for each state, there is an individual nuclear trajectory associated with, and the EOM becomes

To incorporate electronic transitions into the MST method, we invoke the MM model and the classical MM Hamiltonian is given by^2,22

with H_kk and H_kl the nuclear Hamiltonian matrix elements, $x , p$ and $Q , P$ the phase points for the classical electronic and nuclear DOFs, respectively. $n k ≡ 1 2 x k 2 + p k 2 −γ$⁠, and n_kl ≡ x_kx_l + p_kp_l represent the population of the electronic state k and the coherence between state k and l, respectively. $γ= 1 2$ in the quantum representation and can be taken as a parameter in the semiclassical or quasiclassical simulations.³⁸ The electronic dynamics is therefore governed by MM Hamiltonian Eq. (6), i.e.,

The nuclear dynamics in Eq. (5) can be described by the IVR method.^35,36 If we take the nuclear eigenfunctions to be coherent states, i.e., $χ j l = P j l ( t ) Q j l ( t )$⁠, Eq. (5) thus can be rewritten in terms of the IVR,

By using the linear approximation,²⁹ we obtain the following

where $Q ̄ j k l ( t )$ is the average coordinate of the nuclear trajectories on the states j and k. The overlap function in Eq. (8b) could decay quickly as time evolves, then the interstate couplings might be negligible, and nuclear motions reduce to a BO-like dynamics, i.e.,

In this work, we focus on the quasiclassical approximation to the MST approach, i.e., assuming the multi-state nuclear trajectory is classical. The EOM for the lth nuclear trajectory associated with the state j is given by

where F is damping function accounting for the effect due to the separation of the trajectories on different states ΔQ_jkl(t). When the inter-state couplings are small, the individual trajectories are largely independent, which may provide a better description on the multi-state dynamics than the Ehrenfest model. The coupled nuclear-electron dynamics therefore can be described by the ensemble average of the corresponding time-dependent variables [e.g., A(t)] over initial conditions on electronic DOFs and multi-state nuclear trajectories, i.e.,

with N, and n the number of the nuclear and electronic DOFs, respectively; and $ρ x 0 , p 0 , Q 0 , P 0$ the partition function.

The electronic population associated with the multi-state trajectory is determined in the same way as for the SQC approach.³⁸ Specifically a symmetrical binning window function is applied to project the ensemble of multi-state trajectories onto the corresponding electronic state, i.e.,

where Δn = 2γ and h(z) the Heaviside function. The electronic quantum number N taking to be either 0 or 1 corresponds to the unoccupied or the occupied state, respectively. Once the electronic population for the multi-state trajectory is obtained from Eq. (11), dynamical properties depending on electronic population can be evaluated by Eq. (10) using the symmetrical windowing function. For example, time-dependent electronic population of state k is given by the following Monte Carlo average, assuming that the system starts from the state i,

Eq. (12) appears to be analogous to that in the original SQC approach. However the electronic dynamics here is coupled with multi-state trajectories instead of an individual one.

We take the spin-boson system and two Tully’s nonadiabatic models for single and dual avoided crossings as benchmark models³⁸ for testing our MST method. First is the spin-boson model for condensed phase nonadiabatic dynamics determined by the Hamiltonian as follows:

with the bath contribution given by

and

Here F = 100 is the number of bath modes coupled with electronic states, and the frequency and coupling constants ω_kj,  c_kj obey the distribution of the spectra density,

We choose the Ohmic bath given by

and the discretization scheme suggested by Wang et al.⁴⁴ The initial coordinates and momenta of nuclear DOFs are generated in a Monte Carlo scheme from the Wigner transform of the Boltzmann operator. See the references for more details.^38–40 There are four cases for the spin-boson model³⁸ to be considered in this work. They are as follows:

The off-diagonal Hamiltonian matrix elements in the spin-boson model are constants, therefore all inter-state coupling terms in Eq. (9) vanish. If the further approximation is made to assume that the individual nuclear trajectory associated with each state is independent of electronic population, nuclear trajectories are just like those in BO dynamics. We denote MST-BO for this approximate MST treatment and MST-QC for the original quasiclassical MST method in Eq. (9). Note that in MST-BO, electronic population dynamics is still coupled with nuclear DOFs through the nuclear Hamiltonian matrix, as shown in Eq. (7).

The other two benchmark systems are the one dimensional single and dual avoided crossing models used by Tully in the surface hopping treatment for nonadiabatic dynamics.⁶ The Hamiltonian matrix elements for the single avoided crossing problem is given by

where A = 0.01, B = 1.6, C = 0.005, and D = 1. And for the dual avoided crossing system, we take

where A = 0.1, B = 0.28, C = 0.015, D = 0.06, and E = 0.05.

The off-diagonal Hamiltonian matrix elements in the two Tully’s models are nuclear-coordinate dependent, and the inter-state couplings in Eq. (9) are nontrivial. Here we test two different simple choices for the damping function as follows:

and

We denote the MST treatments based on these two choices as MST-QC and MST-QC-linear, respectively.

For interesting readers, the algorithm for the implementation of the MST-QC approach (other versions of the MST approach should be in a similar style) to an n-state nonadiabatic system may be coded in the following steps:

1. Genaerate random numbers to obtain classical electronic variables (x_i, p_i) for each electronic state and classical nuclear phase space variables (Q, P) according to their initial distributions.

2. Propagate classical trajectory dynamics forward with a time step dt for each trajectory associated with the corresponding state according to the EOMs Eq. (7) for electronic DOFs and Eq. (9) for the nuclear DOFs.

3. Determine which electronic state the system is on using the symmetrical windowing function Eq. (11). Note that here even though the total system may be in only one state, all the multi-state individual trajectories are continuously propagated no matter to which state the system belongs.

4. Update the Hamiltonian matrix elements and forces to continue on the next step.

We first examine the spin-boson problem. The initial state of the system is placed on the state 1, and we calculate the time-dependent electronic population difference between state 1 and state 2, i.e.,

Figure 1 displays the D(t) functions obtained by using the MST-QC and MST-BO methods along with the results from the exact quantum calculations (taken from Ref. 38), the SQC approach, and the linearized semiclassical (LSC) IVR method.²⁹ 

Figure 1(a) is for the symmetric unbiased spin-boson problems at high temperature, i.e., case I in Eq. (17), and the corresponding low temperature results for case II are shown in Figure 1(b). The MST-QC approach provides results in excellent agreement with the exact quantum calculations. Even the more approximate MST-BO method describes the coupled nuclear-electron dynamics quite accurately for these two cases. Note that the MST-BO method is not the same as the conventional BO dynamics. In fact, it represents for the opposite limit for which nuclear coherence decays more quickly than electronic coherence, and electronic dynamics in MST-BO is still coupled with nuclear motions which allows for nonadiabatic transitions. The long-lived electronic coherence might be the reason why classical-like descriptions such as SQC and LSC-IVR on the nonadiabatic dynamics in such spin-boson problems work so well.³⁴ 

For the more challenging asymmetric spin-boson problems, as shown in Figures 1(c) and 1(d) for case III and case IV, the MST-QC method predicts the short-time dynamics, oscillations, and high temperature behavior (see Figure 1(d)) quite faithfully, but a slow decay rate at low temperature (see Figure 1(c)). For case III, the SQC approach is the best and the MST-QC method is as good as the LSC-IVR method except for stronger coherence. The MST-BO treatment reproduces the decay rate poorly, probably due to the insufficient inclusion of the “quantum backreaction.”² Choosing a smaller value (0.366) for the parameter γ results in stronger coherence, and the results are presented in Figure 2. The effect of choosing different “optimized” gamma values [0.366 vs 0.5, see discussions in Ref. 38] appears to be small for this model and the following other models in this work. However the check on this effect is suggested to be done for other general applications.

Now we test the MST method on Tully’s models for nonadiabatic dynamics. Figure 3 displays the transmission probability of the nuclear wavepacket passing through the single avoided crossing of two electronic states. Taking γ = 0.5, the MST-QC approach provides the transmission probabilities in good agreement with the quantum exact results in a wide range of initial nuclear momentum except for the very low momentum region. The MST-BO results are a little worse than the MST-QC calculations. For γ = 0.366, both MST-QC and MST-BO predicts a higher T_2←1 for the intermediate nuclear momentum but a little lower T_2←1 for large nuclear momentum. The overall agreement with the exact results does not change much with the adjustment of γ. It seems that the MST-QC works a little bit better for γ = 0.5, whereas the MST-BO appears better for γ = 0.366.

We then check the effect of using a different damping function given by Eq. (20b) in the EOM of nuclear dynamics. The results calculated from the MST-QC and MST-BO methods along with their corresponding linear treatments for the damping function are shown in Figure 4. Introducing an additional linear term with respect to the separation between multi-state nuclear coordinates into the damping function produces nearly the same transmission probabilities, especially for large nuclear momentum. Many other choices for damping function could be taken to represent for inter-state nuclear coherence more accurately, which might improve the MST predictions in the low nuclear momentum region. However, at least in the current applications (also refer to the dual avoided crossing model in Figure 6), the zero order classical approximation given by Eq. (20a) seems adequate.

The dual avoided crossing problem is more complicated and exhibits the Stuckelburg oscillation in transmission probability. The MST-QC and MST-BO results are presented in Figure 5 along with the exact quantum calculations. In this case, the MST-QC method largely reproduces the oscillation; however, the predicted magnitude and the phase of the oscillation are not quantitatively good. The MST-BO overestimates the magnitude of the oscillation and appears worse than the MST-QC treatment. Changing the value of parameter γ from 0.366 to 0.5 does not improve the results substantially.

Figure 6 displays the MST-QC and MST-BO results in comparison with those obtained by taking a different damping function given by Eq. (20b). Again the effect is negligible in this case. The linear version of MST-QC method produces almost identical transmission probabilities in the entire energy region. For the MST-BO method, only small change in low energy region is noticed.

It is worth mentioning that the SQC approach has been demonstrated to be quite accurate and efficient in a variety of nonadiabatic applications^37–42 including the benchmark models considered in this work.³⁸ However the nuclear dynamics in the SQC is still Ehrenfest-like, that is, to say the nuclear dynamics of any individual trajectory is governed by an average of multi-state forces. Alternatively the MST method propagates a multi-state trajectory simultaneously, and for each trajectory associated with a particular state, a multi-state force is applied (see Eqs. (5) and (9)). Unlike the averaged force in the SQC approach or other Ehrenfest-like dynamics, inter-state nuclear coherence is explicitly built-in for the multi-state force in the MST description. Therefore for the nonadiabatic dynamics involving many states, advantages of the MST treatment in being able to shut off nuclear coherence for the insignificant DOFs or completely for all nuclear DOFs may achieve neat physical pictures and excellent numerical stability. Presumably the multi-state trajectory in MST could represent nuclear dynamics for each state more faithfully than the individual trajectory in Ehrenfest-like dynamics or the state-switched trajectory in surface-hopping dynamics. Keep in mind that here only the ensemble average of the multi-state trajectory makes physical senses when comparing with the classical pictures.

Moreover when the further approximation is made to turn off the electronic dependence in the nuclear dynamics for the individual nuclear trajectory associated with each state, the resulting MST-BO method is very simple and easy to be implemented. This is in particular important for treating nonadiabatic dynamics in realistic molecular systems involving many states. The implementation of the MST-QC and MST-BO methods to electron coupled lithium diffusion in a commercialized lithium ion battery material LiFePO₄ has been done,⁴³ and the work will be published elsewhere.

By construction, the MST method we suggest here normally requires multiple nuclear trajectories to be simulated, therefore the computational cost for the nuclear dynamics is scaled with the number of states. However the overall performance is still encouraging, since the dynamics is classical. For the current applications presented in the paper, 960 000 trajectories (96 electronic times 10 000 nuclear configurations) are used for Tully’s models, and 96 000 trajectories (96 electronic times 1000 nuclear configurations) are for the spin boson problem. However using ten times fewer trajectories can still provide very good results. Figures 7 and 8 provide the convergence check on the model systems considered in this work for the MST-QC method. The results for the MST-BO method are similar (not shown). For the spin boson problem, the result obtained by using only 10 000 trajectories (10 electronic times 1000 nuclear configurations) converges very well (case III, Figure 7). Only small deviations come out for using the nuclear configurations as few as 100 (96*100). The convergence check on other cases in Figure 1 gives similar results (not shown). For the Tully’s models, reducing the number of electronic configurations (from 96 to 10) or nuclear configurations (from 10 × 10³ to 1 × 10³) results in no noticeable changes in the probability calculations. Even for the results by using only 10 000 trajectories, the largest deviation is less than a few percents.

A MST method is proposed to provide a trajectory representation for the coupled nuclear-electron dynamics. In the MST method, the nuclear motions evolve in terms of a multi-state trajectory with each state associated with an individual trajectory, and the coupled electronic dynamics is described by the MM model. The QC and BO approximations of the MST approach have been tested on the spin-boson problem and two Tully’s models for nonadiabatic dynamics. The MST-QC treatment provides results in fairly good agreement with exact quantum calculations for the benchmark models. For the symmetric spin-boson problem and the single avoided crossing model, even the further approximate MST-BO works reasonably well. Further improvements on the MST method might be made by a better description on the interstate nuclear coherence, for example, to introduce semiclassical individual trajectory for the multi-state trajectory, and to optimize the damping function, etc. The development of the ab initio MST approach to perform nonadiabatic quantum dynamics simulations on-the-fly seems also promising. The MST method and its BO approximation in particular have the potential to be implemented to nonadiabatic dynamics in realistic complex molecular systems. Also the MST method is not only for electronically nonadibatic dynamics but also should be able to be applied to any general nonadiabatic systems such as vibrationally nonadiabatic problems.

We acknowledge the support from Peking University Shenzhen Graduate School, Shenzhen Science and Technology Innovation Council (Nos. JCYJ20120829170028565 and ZDSY20130331145131323), National Science Foundation of China (No. 51471005), and National Supercomputing Center in Shenzhen (Shenzhen Cloud Computing Center).

The authors declare no competing financial interest.