Nonadiabatic Dynamics Algorithms with Only Potential Energies and Gradients: Curvature-Driven Coherent Switching with Decay of Mixing and Curvature-Driven Trajectory Surface Hopping

Direct dynamics by mixed quantum–classical nonadiabatic methods is an important tool for understanding processes involving multiple electronic states. Very often, the computational bottleneck of such direct simulation comes from electronic structure theory. For example, at every time step of a trajectory, nonadiabatic dynamics requires potential energy surfaces, their gradients, and the matrix elements coupling the surfaces. The need for the couplings can be alleviated by employing the time derivatives of the wave functions, which can be evaluated from overlaps of electronic wave functions at successive timesteps. However, evaluation of overlap integrals is still expensive for large systems. In addition, for electronic structure methods for which the wave functions or the coupling matrix elements are not available, nonadiabatic dynamics algorithms become inapplicable. In this work, building on recent work by Baeck and An, we propose new nonadiabatic dynamics algorithms that only require adiabatic potential energies and their gradients. The new methods are named curvaturedriven coherent switching with decay of mixing (κCSDM) and curvature-driven trajectory surface hopping (κTSH). We show how powerful these new methods are in terms of computer time and good agreement with methods employing nonadiabatic coupling vectors computed in conventional ways. The lowering of the computational cost will allow longer nonadiabatic trajectories and greater ensemble averaging to be affordable, and the ability to calculate the dynamics without electronic structure coupling matrix elements extends the dynamics capability to new classes of electronic structure methods.


Introduction
Simulating the dynamics of processes involving electronically excited molecules requires potential energy surfaces, their gradients, the electronic matrix elements controlling their coupling, and a nonadiabatic dynamics algorithm. 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19 Recent advances on these topics have enabled quantitative accurate simulations for medium-size molecules and qualitatively accurate simulations for nano-sized molecules or clusters. 20,21,22,23,24,25,26,27,28,29,30,31 The electronic structure must be treated quantum mechanically, and for simulating all but the smallest systems, one uses a semiclassical treatment of nuclear motion. The present article is concerned with direct dynamics, where, instead of requiring parametrized analytic functions for the energies and couplings, all required energies, forces, and couplings for each geometry that is important for evaluating dynamical properties are obtained directly from electronic structure calculations when they are needed in the dynamics calculation.
Electronic structure calculations directly yield electronically adiabatic wave functions, ! , and the coupling in the semiclassical calculations in the adiabatic representation is provided by the nuclear momentum operator, which involves the matrix elements of the gradient with respect to nuclear coordinates R: 1,2,32 (1) where denotes an electronic matrix element. The matrix element in eq 1 is usually called the nonadiabatic coupling (NAC). It is a 3N-dimensional vector where N is the number of atoms. The present article is concerned with an efficient approximation to the NAC.
In semiclassical methods, the electronic wave functions may also be considered to be functions of time (t) because (2) where r denotes the electronic coordinates, and R(t) is the trajectory of nuclear motion. A popular semiclassical method is trajectory surface hopping (TSH). 1,2,9,33,34,35,36,37,38 In TSH one only requires the NAC in the direction of the current velocity, ̇, because the semiclassical equations only require the time derivative coupling (TDC): (3) IJ Employing a direct calculation of the TDC is more efficient than calculating the NAC, and it can be shown to yield more numerically accurate integration of the TSH electronic equations of motion than direct use of the NAC, especially for situations of trivial crossings. 39 where and are adiabatic potential energies, and q is a one-dimensional nuclear coordinate.
The CSDM algorithm is based on SE, but in CSDM the SE algorithm is augmented by non-

Theory
Since eq 4 is the one-dimensional derivative with respect to a coordinate in the direction of motion, we generalize it to the following one-dimensional derivative with respect to time: where is the local gap between adiabatic potential surfaces, and we use κ as a prefix and as a superscript to denote approximation based on the curvature of the gaps. Since the NAC is skew-Hermitian, we also have Employing eqs 5 and 6 for the TDC is the central feature of κSE and κCSDM. This approximation can also be employed in TSH, yielding κTSH. We have implemented κSE, κCSDM, and κTSH in a new version of SHARC-MN, which is available at our cost at our software distribution site. 72 Next we give (i) the details of the how we calculate the TDC numerically in SHARC-MN and (ii) the remaining details necessary to completely specify κSE and κCSDM, which -as stated above -require an effective NAC as well as the TDC.
κTDC. The curvature-approximated TDC is given in general by eqs 5 and 6. However, there are some details that need to be clarified.
First, when the radicand is negative, we simply set it to zero. This is justified as follows. The reason why a formula like eq 4 or q 5 works is best appreciated by considering the nature of two potential curves near a locally avoided crossing. In a typical case the lower curve is concave, and the upper curve is convex, so the radicand is positive. Therefore, a situation where the radicand is negative is expected to occur only far from conical intersections where the NAC is small, so this should not cause a large error.
Second, we note that can be re-written as, In the practical SHARC-MN implementation, this is computed from backward finite differences as follows: where (9) SHARC employs velocity-Verlet algorithm, 73 in which the velocity at time step is evaluated after the nuclear force at is computed.
For both SE and CSDM, the nuclear forces are not decidable until one has determined the electronic coefficients. And therefore, , which is needed when eq 9 is substituted into eq 8, is not available for evaluating when one uses the velocity-Verlet algorithm.
Hence, in the current implementation, is approximated by a forward propagation, where is acceleration vector at time t.
For TSH, when uses the velocity-Verlet algorithm, one does have when evaluating the TDC. However, when a hop happens, the velocity needs to be adjusted to enforce energy conservation. Therefore, to make a consistent approximation over the whole trajectory, we employ eq 10 to evaluate for κTSH as well.
where is the coefficient of adiabatic electronic state . Inserting eq 11 into the time dependent electronic Schrödinger equation, using eqs 5 and 6 for the TDC, and performing some straightforward calculus gives the electronic equation of motion (EOM) for κSE: The nuclear EOM of the regular SE method is (13) where P is the momentum conjugate to R, and is an element of the electronic density matrix. For the tSE, method, we re-derived the nuclear EOM such that it depends on an effective NAC: (15) where is the effective NAC defined by the difference gradient vector , the velocity vector Ṙ, and the TDC: We justified calling an effective NAC because it satisfies (19) We now replace in eq 18 by , such that eqs 17-19 yield a new effective NAC, which will be called the curvature-approximated effective NAC and denoted by .
As before, 75 it is necessary to project out the rotational and translational components of the effective NAC in order to conserve the total angular momentum and the position of the center of mass. This yields 75 (20) where 1 and Q are the identity matrix and a projection operator. 75 The operation in eq 20 removes unphysical translational and rotational components of . Therefore the nuclear EOM of κSE becomes: Equations 12 and 21 define the electronic and nuclear EOM for κSE.
κCSDM. Adding decoherence to κSE in the same that it was added to tSE to get tCSDM 53 yields the following electronic and nuclear EOMs for κCSDM:

( )( )
Re IJ   I  J  IJ  I  II  IJ  I  IJ I where and denote the electronic and nuclear equation of motion from κSE, as given in eqs 12 and 21, and where the additional terms are decay-of-mixing (DM) terms, which are explained in detail elsewhere. 48,53 In tCSDM decoherence contribution to the change in the nuclear momentum is 53 (24) in which is the decoherence time, K is the pointer state, and is the decoherence vector for states I and K given by (25) where is the internal vibrational momentum, and a0 ≡ 1 bohr.
In CSDM, tCSDM, and κCSDM, we propagate two density matrices; the true density matrix is propagated by equations 22 and 14, and another density matrix, called the coherent density matrix, is propagated coherently (i.e., by eqs 12 and 14) and is used to control the switching of the pointer state. Although the coherent density matrix is in general different from true density along the trajectory, it is re-initialized to true density for every complete passage of a strong interaction region.
Beyond the changes already discussed for κSE, two additional changes are required to obtain κCSDM: In tCSDM, the complete passage of a strong interaction region is defined by local minima of (26) In κCSDM, we replace this by (27) The decoherence direction of eq 25 is replaced by

Computational details
We implemented κSE, κCSDM, and κTSH in SHARC-MN, 72 and we tested κCSDM against the previously published CSDM and tCSDM methods. We used nonadiabatic dynamics of the ethylene molecule as our test example. The Molpro software package 76 was interfaced to SHARC-MN to perform direct dynamics.
Ground-state geometry optimization of ethylene together with vibrational analysis were performed with MP2/6-31G**, 77,78 and the initial geometries and velocities of the trajectories were randomly sampled from the ground-vibrational-state Wigner distribution. Before propagating the EOMs, the molecule was raised vertically to the first excited state (S1).
We performed tests of κCSDM and κTSH. We do not test κSE because we use SE and κSE simply as base methods upon which to build CSDM and κCSDM -not preferred choices for practical simulations.
The potential energy surfaces and gradients for ethylene direct dynamics were calculated by state-averaged complete-active-space self-consistent-field theory (SA-CASSCF) 79,80 by averaging over 3 states with a (2,2) activated space (two active electrons in two active orbitals) and the 6-31G** basis set.
We ran 300 trajectories for κCSDM and 150 trajectories each for CSDM and tCSDM. We We also ran trajectory surface hopping with the energy-based decoherence correction 81 (TSH-EDC). We ran 200 TSH-EDC trajectories -100 using the time derivative calculated as an overlap and 100 using the scalar product of the NAC and the velocity. With converged ensemble averaging and converged step sizes and if one sued the same velocity rescaling at hops, one would get the same result by these two methods so in that sense it is legitimate to combine their results; however, we ran these trajectories with both algorithms so we could compare the timing.
When we need to distinguish thee two algorithms, we label them as TSH-EDC(TDC) and TSH-EDC(NAC).
Finally, we ran 100 trajectories for κTSH-EDC. The decoherence correction is the same as in TSH-EDC. The velocity component in the direction of the gradient difference is rescaled to adjust the kinetic energy after a successful surface hop for κTSH-EDC and TSH-EDC(TDC), but for TSH-EDC(NAC), the direction of the projected NAC is used to adjust kinetic energy after a successful surface hop.

Results
The present results for κCSDM are compared to our previous results 53 for CSDM and tCSDM. To assess the excited-state lifetime, we employ ensemble-averaged population analysis using the following definition of the probability of being in state I: (29) where is the index of trajectory, is the total number of trajectories, and the probability is where "#$% ! is the number of trajectories for which the pointer state is I when the trajectory is terminated. This population is denoted as the pointer-state population. The pointer-state population for κCSDM, tCSDM, and CSDM is shown in Fig. S1 in Supporting Information.
The ensemble-averaged potential energies for S0, S1, S2 and the effective PES (eq 8 of ref 48) are shown in Fig 2. The effective PES (denoted as Veff) is gradually changing from V1 to V0 around 50 fs which is consistent with the time region where population transfers from S1 to S0.
The half-life (τ1/2) obtained from the simulations for decay of the S1 excited state are 60, 52, and 51 fs for κCSDM, tCSDM, CSDM, respectively, which shows very good agreement among the three methods; in fact, the difference is of the same order as the statistical uncertainty in the halflife due to the finite amount of Monte Carlo ensemble averaging with 150-300 trajectories). To assess the accuracy of curvature-approximated TDC and the curvature-approximated effective NAC, we randomly picked three κCSDM trajectories and performed NAC calculation at each geometry along the trajectory. The norm of ab initio NAC and projected curvatureapproximated effective NAC for S0 and S1 are shown in Fig. 3. For comparison, this figure also shows V0 and V1 along the κCSDM trajectory. It is clear that and spike at the same position, which is where V0 and V1 get close. However, we also observe a trend that seems to be relatively smaller than , and that may explain why population transfer from S1 to S0 seems to be slightly faster in tCSDM and CSDM than in κCSDM; however, it is not clear if the difference is statistically significant. The bottom line is that the curvature-driven method is in good agreement with tCSDM and CSDM.
The computational cost ratio for κCSDM, tCSDM, and CSDM is roughly 1:3.5:4.5. These ratios show a significant speedup for κCSDM, which occurs because no overlap integrals or nonadiabatic coupling vectors are computed from electronic structure software.
Next, we discuss κTSH. The S0, S1, and This shows that the speedup of κTSH-EDC relative to TSH-EDC(NAC) is more than a factor of two even for this small molecule (as the molecule simulated gets larger, this timing savings is expected to become more significant; it could even be dramatic). However, κTSH is only 30% faster than the more accurate κCSDM method.
Furthermore, these timings show that κTDC drops the cost of the mean-field, dual-density-matrix κCSDM algorithm below that of TSH-EDC(NAC).

Concluding Remarks
In this work, we proposed a curvature-approximated time-derivative coupling formula, and we showed how to use this formula to obtain curvature-driven approximations to trajectory    in red and green respectively. We also show the adiabatic potentials V0 (yellow) and V1 (blue) for comparison. Note that states S0 and S1 are numbered 1 and 2 in the row and column indices of the matrix elements.