A multiple time step algorithm for trajectory surface hopping simulations

A multiple time step (MTS) algorithm for trajectory surface hopping molecular dynamics has been developed, implemented, and tested. The MTS scheme is an extension of the ab initio implementation for Born–Oppenheimer molecular dynamics presented in [J. Chem. Theory Comput. 14, 2834 (2018)]. In particular, the MTS algorithm has been modified to enable the simulation of non-adiabatic processes with the trajectory surface hopping (TSH) method and Tully’s fewest switches algorithm. The specificities of the implementation lie in the combination of Landau–Zener and Tully’s transition probabilities during the inner MTS time steps. The new MTS-TSH method is applied successfully to the photorelaxation of protonated formaldimine, showing that the important characteristics of the process are recovered by the MTS algorithm. A computational speed-up between 1.5 and 3 has been obtained compared to standard TSH simulations which is close to the ideal values that could be obtained with the computational settings considered.


I. INTRODUCTION
Non-adiabatic phenomena such as photo-physical or photo-chemical processes are characterized by a failure of the Born-Oppenheimer (BO) approximation commonly invoked to describe molecular systems. The BO approximation, or the closely related adiabatic approximation, decouples the description of the nuclei and electrons. In non-adiabatic processes this coupling becomes important and the BO approximation breaks down.
Many different methods have been developed in order to describe non-adiabatic processes, ranging from fully quantum and formally exact models to mixed quantum/classical or semiclassical approaches (for a review see Refs. 1 and 2). Each method has its pros and cons, but in most cases it boils down to a compromise between computational efficiency and accuracy.
In this work, we focus on one of the most popular methods, the trajectory surface hopping (TSH) approach.
In the TSH method (summarized in section II A), the evolution of the system is represented by a swarm of independent classical nuclear trajectories, which can hop from one electronic state to another in a stochastic way. The forces acting on the nuclei are calculated on-the-fly along each trajectory and transitions between electronic levels are considered simultaneously. In this way, the TSH method is thus able to describe non-adiabatic phenomena such as photo-chemical and photo-physical processes.
The TSH approach belongs to the mixed quantum/classical class of methods and is one of the computationally most expedient way to include non-adiabatic effects. Nonetheless, TSH simulations require the evaluation of the nuclear forces from first-principles simulations, i.e., by solving the time-independent electronic Schrödinger equation, and those forces have In recent years, several attempts to extend the application range of non-adiabatic MD have been proposed. [7][8][9] The present article is a further contribution to reduce the computational cost of non-adiabatic dynamics.
The strategy investigated in this work, is to reduce the computational requirements of the TSH method by relying on a multiple time step (MTS) algorithm. MTS techniques have first been introduced by Tuckerman et al. in the context of classical MD. 10 The MTS scheme relies on a decomposition of the atomic or nuclear forces into different components with different characteristic time scales. This decomposition enables to calculate the slow components of the forces less frequently than the fast one, while maintaining a fully timereversible symplectic propagation. If the computational cost of the slow components is significant, large computational speed-ups can thus be obtained.
This article is organized as follows. After introducing an MTS algorithm for TSH simulations in section II B, the new method is applied to the photorelaxation of protonated formaldimine (section III). The MTS-TSH algorithm is compared to standard TSH simulations both in terms of accuracy and computational cost. Finally, some concluding remarks and perspectives are given in section IV.

II. THEORY
In this first section, we briefly review the TSH formalism for non-adiabatic MD with particular emphasis on the version implemented in the CPMD plane-wave package. 11

A. Trajectory surface hopping in CPMD
The TSH method can be seen as an attempt to introduce coupling between the electronic and nuclear degrees of freedom in Born-Oppenheimer molecular dynamics (BOMD). 12 Standard BOMD consists in a description of the nuclear coordinates based on classical mechanics, where the index α denotes a given nucleus of mass M α and classical coordinates R α and the two dots on top of the coordinates denote a second-order derivative with respect to time. When no index is specified, R stands for all nuclear coordinates. The forces acting on the nuclei, F α (R), are usually calculated on-the-fly, from ab initio electronic structure calculations at fixed nuclear geometries, where the diagonal matrix elements of the molecular electronic Hamiltonian,Ĥ (under the BO approximation) are given by, In Eq. (3) we have introduced the electronic wavefunction, ψ L , for an arbitrary adiabatic state, L, which depends on the electronic coordinates r. The parametric dependence of the electronic wavefunction on the nuclear geometry, R, is also given. Generally, a single BOMD trajectory will thus propagate the classical nuclear coordinates on a single PES corresponding to a specific adiabatic electronic state.
However, when more than one electronic state is important to describe the dynamics of a system (for example in the description of photo-physical phenomena) it is important to go beyond the BO approximation and consider non-adiabatic algorithms.

Tully's fewest switches method
One of the most popular approaches used to describe such phenomena is the TSH method, in particular when combined with Tully's fewest switches algorithm. 2,13,14 In TSH a given trajectory can hop from one electronic state to another in a stochastic way, depending on the probability of the transition to occur. In order to get the transition probabilities one generally has to solve a time-dependent equation for the electrons of the system, Where the time-dependent coefficients C J (t) comes from an expansion of the time-dependent electronic wavefunction as a linear combination of time-independent adiabatic states. N states is the total number of electronic states considered, ω J is the excitation energy for state J and σ JK is the non-adiabatic coupling (NAC) term, In CPMD, the NAC terms are computed by finite differences and using a CIS representation of the excited states. [15][16][17] Finally, in the fewest switches (FS) scheme, the probability of transition from adiabatic states J to K is evaluated as, where δt is the classical time step used to integrate eq. (1), R(z) denotes the real part of z and negative probabilities are set to zero. The decision to hop from state J to state K is then taken by generating a random number r ∈ [0, 1] and evaluating the following condition,

Landau-Zener transition probabilities
As a simpler alternative to the solution of eq. (4) for the calculations of the transition probabilities in eq. (6), it is possible to obtain approximate transition probabilities from Landau-Zener-Stückelberg (LZ) theory for non-adiabatic transitions. 18,19 In the CPMD package, such probabilities are computed directly from the knowledge of the energy of the adiabatic electronic states as, where ∆E adia JK is the gap between adiabatic states J and K directly obtained as a byproduct of DFT and TDDFT calculations. 18 In a TSH simulation relying on LZ theory, a hop from an electronic state to another is considered based on the following condition, where r is again a random number chosen between zero and one.
For more details about the implementation of TSH in the CPMD package, see Refs.
where L is the Liouville operator given by, and F j is a single component of the nuclear forces (the index j is a collective index for Cartesian coordinates and nuclei). By assuming a time scale separation of the forces into fast (F fast ), and slow (F slow ) components, it is possible to rewrite the Liouville operator as, iL slow Applying a Trotter factorization on the classical propagator and discarding terms of third order and higher in t, we can define a discrete time propagator, which can be translated into an MTS algorithm, 10,22 with G fast (∆t/N ) = G fast (δt) = e iL fast p (δt/2) e iLxδt e iL fast p (δt/2) .
In eqs. (16) and (17) we have introduced two finite time steps; ∆t, which reflects the time scale of the slow forces (F slow ), and δt = ∆t/N , which is adapted to the fast forces (F fast ).
It should now be apparent that such an MTS algorithm can lead to computational savings, if the slow forces that have to be calculated less frequently are computationally more demanding.
Algorithm 1 represents a pseudo-code that can be obtained by applying each term of the propagator in eq. (16) (one by one from the right to the left) onto an initial phase space element Γ(0) ≡ {x, p}. In the context of first-principles BOMD, the separation of the forces into fast and slow components is not straightforward. In the CPMD package we have chosen to use different levels of electronic structure theory to decompose the forces. Typically, a "low" level density functional (e.g. GGA) is used to calculate the fast components of the nuclear forces, while the slow components are obtained as the difference between the forces obtained with a "higher" level functional (e.g. hybrid) and the "low" level forces, The physical motivation for this separation comes from the fact that the chosen electronic structure levels differ only in their treatment of correlation (or exchange and correlation in the case of DFT). Since those contributions correspond to relatively weak interactions in terms of energy (with no explicit dependence on nuclear positions) they can be expected to represent relatively weak/slowly varying force contributions.
The implementation of the MTS algorithm in CPMD makes use of this separation as well as a slightly different but completely equivalent layout of the code. 21 This implementation presented in pseudo-code in algorithm 2, makes the MTS algorithm look like a velocity Verlet algorithm with effective forces, F eff , that are time step dependent.
Even though, the force decomposition in terms of high and low electronic structure levels is done ad hoc, this kind of separation has already proven useful 21,23 and the benefits in terms of computational cost are evident.

MTS algorithm for trajectory surface hopping dynamics
When using the MTS implementation described in algorithm 2 in combination with a TSH algorithm, one has to decide how to hop from one electronic state to another. In order to get trajectories of high accuracy, it would be beneficial to consider electronic transitions based on Tully's FS criterion in eq. (7) calculated with the MTS high level functional. In the following, this type of calculation (with a velocity Verlet algorithm) will actually be used as a reference. However, when using an MTS algorithm, if the outer time step ∆t becomes large, some parts of the PES where the transition probabilities are high might be treated only by the low level functional and the transitions would be missed at the high level.
In this work, we propose another strategy that can potentially solve that problem. This strategy consists in evaluating the transitions probabilities with the LZ formula in eq. (8) 1: Initialize positions, velocities, high and low level forces: x, v, F high , F low 2: Get effective force: Get both high and low level forces, F high , F low 8: Get effective force: else (inner step) 10: Effective forces are set to the low level forces: F eff ← F low 11: end if (outer/inner steps) 12: ALG. 2: MTS-BOMD algorithm as implemented in the CPMD package. This algorithm can be obtained straightforwardly from algorithm 1 by using the partitioning of the forces in eqs. (18) and (19) and reshuffling a few steps.
during the low level steps, while for high level steps, the electronic transition are evaluated according to the FS criterion in eq. (7) based on the high level quantities. This is not yet completely satisfactory since it does not guarantee that a transition detected with the LZ probabilities during a low level step would have also been detected by the FS criterion calculated with the high level functional. To further improve on that issue, we suggest that if a transition is detected at the low level (using LZ theory), a high level calculation is triggered to confirm the transition using Tully's FS criterion. This strategy is described in algorithm 3 and tested in the remaining sections.
Finally, it is important to note that in the case of high level calculations triggered by the low level LZ criterion, the random number used in the high level FS criterion in eq. (7), should be the same as in the low level LZ criterion in eq. (9). We also underline that the discrete time step δt in eq. (6) always correspond to the inner time step in the MTS scheme.
This can be rationalized by realizing that at each inner time step, an electronic transition can occur if it is detected at the low level and confirmed at the high level, such that when eq. (6) is invoked, it is only to check for a transition in the last δt time window.

III. RESULTS AND DISCUSSION
In this section, the protonated formaldimine (CH 2 NH + 2 , denoted as system I) is used as a simple yet interesting example to investigate the capabilities of the new MTS-TSH method presented in section II B. In particular, we will investigate the possibility to use the MTS-TSH algorithm as a more efficient alternative to the standard FS-TSH algorithm. The physical process under investigation in this section is the photorelaxation of system I. This process is a typical example of photo-dynamics and is thus very handy to test new models for non-adiabatic molecular dynamics. 9,15,[24][25][26] All the calculations presented in this section have been performed with a local version of the CPMD plane-wave package. 11 For the reference TSH simulations, the nuclear forces and NACs are computed with the PBE0 hybrid functional. 27 The same functional is thus if FS criterion in eq. (7) is met for any state K = J then if FS criterion in eq. (7) is met for any state K = J then  In Fig. 1 we have represented the evolution of the transition probabilities (P S 2 →S 1 ) for the three different levels. From Fig. 1, it is clear that neither the PBE-LZ probabilities nor the PBE-FS probabilities represent a completely reliable approximation to the reference PBE0-FS transition probabilities. However, a relatively good correlation is observed. In particular, whenever the reference PBE0-FS transition probabilities are large, the PBE-LZ probabilities are also relatively large. For some reason this is less obvious for the PBE-FS curve.
The algorithm developed in section II B relies on the low level calculations to detect potential electronic transitions. Whenever a low level transition is detected, it is double checked with a high level FS calculation such that the low level probabilities do not have to be quantitatively accurate. The transition probabilities in Fig. 1 indicates that using LZ theory at the low level for the calculation of the transition probabilities is enough. In other words, the LZ probabilities can be used as a proxy during the low level steps of the MTS-TSH method.
We note that the objective of the investigation performed in this section is to support the design of the algorithm presented in section II B and that further tests could be performed to draw more general conclusions. Nonetheless, since the usage of LZ probabilities at the low level is only used to trigger high level calculations, we believe that a strong empirical support is not required at this stage.

B. Comparing single trajectories via deterministic surface hopping
As we have seen in section II A, non-adiabatic dynamics performed with a TSH algorithm are stochastic by nature due to the randomness used in the hopping procedure. This stochastic behaviour makes it difficult to compare individual trajectories obtained with a TSH algorithm. To properly compare different non-adiabatic models, one needs to look at a statistical ensemble of trajectories. Before we present such results in section III C, we first consider in this section a deterministic version of TSH in which the random number r used in eqs. (7) and (9) has been fixed arbitrarily to r = 0.3. We note that, the only MTS algorithm tested here and in the next section is the one presented in section II B, i.e., low level (PBE) LZ transition probabilities are used during the inner steps to trigger a high level (PBE0) calculation which confirms or not the electronic transition.

Quality assessment
Six different runs have been produced, all starting in the second excited state and with the same nuclear configuration. For simplicity, the nuclear velocities are initialized to zero.

Efficiency assessment
Let us call t FS = t high the average CPU time per step spent in a standard FS-TSH simulation with a "high" level functional. This time takes into account the SCF optimization, the solution of the TDDFT equations as well as the calculation of the FS probabilities. In the following we use t FS as a reference CPU time. In order to provide fair comparisons, we also need to consider the average CPU time per step from a LZ-TSH simulation with a "low" level functional t LZ = t low .
The expected or ideal averaged CPU time per step in the MTS-TSH algorithm can then be calculated as, where N is the MTS factor and, t high and t low are CPU timings coming from standard (non MTS) simulations. X denotes the average frequency of triggered high level steps which is a quantity difficult to predict. To obtain the ideal MTS-TSH speed-up we set X = 0, and get, In the limit of a negligible cost of the low level steps (compared to the high level ones) we get, In practice, several things can impact the ideal and limit speed-ups. The most obvious one being the number of triggered high level steps. It is easy to realize that for large values of N , the number of triggered high level steps will tend towards a system dependent number, in most cases larger than zero. Such that one cannot achieve arbitrarily large speed-ups simply by increasing N . From a more practical point of view, the physics of the system (vibrational frequencies) will be the first parameter to consider as a limitation for the value of N and thus for the speed-up (too large values of N would lead to artifacts such as resonances 30 ). A less straightforward impact on the effective (or real) speed-up is given by the overhead due to the convergence of the high level (TD)DFT parameters. Indeed, since with increasing values of N , larger nuclear displacements occur between high level steps, the initial guess for the electronic structure calculation becomes less appropriate, often resulting in more (e.g. SCF) iterations and thus higher computational requirements.
In Fig. 5, we have represented the ideal and limit MTS speed-ups as calculated from eqs. (21) and (22) From Fig. 5, we can see that both the real and ideal speed-ups are quite far from the limit speed-ups (up to a factor 3 of difference). This is simply a consequence of the fact that the condition, t high t low is not satisfied here. For δt = 10 a.u. t high 5.1 · t low , while for δt = 15 a.u. t high 4.8 · t low . It means that one way of improving the efficiency of MTS techniques is to consider cheaper "low" level models.
The real and ideal speed-ups are much closer from each-other and with an MTS factor between 3 and 6 it seems that reliable results could be obtained with a speed-up factor ranging from 1.  In Fig. 6, the collective evolution of the state populations along the dynamics are represented. The main characteristics of the photorelaxation process are well reproduced by all three MTS-TSH runs. In particular, the population is transferred from S 2 to S 1 in the first 10 to 20 fs and then slowly decays into the ground state. We note that the tail of the population of S 2 seem to become larger with the MTS factor, which shows the limit of the MTS scheme. This is reflected in the average lifetime of S 2 , for which we get 8 shows that different decay mechanisms are possible, hopping directly from S 4 to S 2 and then to S 1 or hopping first to S 3 and then directly to S 1 . Those mechanisms are rare due to the fact that the initial population of S 4 is much lower than the population of S 2 , but they are also part of the MTS-TSH swarm of trajectories, which indicates that the new MTS approach is reliable.
Following the work of Westermayr et al. in Ref. 9, we have analyzed the geometries at which the S 2 to S 1 and S 1 to S 0 transitions occur. Fig. 7 represents the values of relevant geometrical parameters at the hopping geometries for the first and second transitions, respectively. From Fig. 7, we can see that the hopping geometries from the MTS-TSH runs spread basically over the same region as the hopping geometries from the reference TSH simulations.
As the MTS factor is increased, the average number of standard high level steps per trajectory will decrease, which should decrease the computational cost of the MTS-TSH   3.12 for the MTS factors 4, 6, and 8, respectively, as can be seen from Fig. 8.
Overall, this statistical investigation indicates that the MTS-TSH algorithm introduced in section II B allows to reproduce results from standard TSH simulations with a significant speed-up.

IV. CONCLUSIONS AND OUTLOOK
We have presented a new algorithm for non-adiabatic molecular dynamics simulations This is achieved by pre-evaluating the transition probabilities during inner steps with a low-level LZ criterion. If a transition is detected, a high level calculation is triggered, to confirm (or not) the electronic transition.
This new MTS-TSH algorithm has been tested successfully on the photorelaxation of protonated formaldimine. We have shown that the MTS-TSH method is able to recover the correct state population along the reaction path as well as the correct geometries at