Adaptive time steps in trajectory surface hopping simulations

Trajectory surface hopping (TSH) simulations are often performed in combination with active-space multi-reference configuration interaction (MRCI) treatments. Technical problems may arise in such simulations if active and inactive orbitals strongly mix and switch in some particular regions. We propose to use adaptive time steps when such regions are encountered in TSH simulations. For this purpose, we present a computational protocol that is easy to implement and increases the computational effort only in the critical regions. We test this procedure through TSH simulations of a GFP chromophore model (OHBI) and a light-driven rotary molecular motor (F-NAIBP) on semiempirical MRCI potential energy surfaces, by comparing the results from simulations with adaptive time steps to analogous ones with constant time steps. For both test molecules, the number of successful trajectories without technical failures rises significantly, from 53% to 95% for OHBI and from 25% to 96% for F-NAIBP. The computed excited-state lifetime remains essentially the same for OHBI and increases somewhat for F-NAIBP, and there is almost no change in the computed quantum efficiency for internal rotation in F-NAIBP. We recommend the general use of adaptive time steps in TSH simulations with active-space CI methods because this will help to avoid technical problems, increase the overall efficiency and robustness of the simulations, and allow for a more complete sampling.


I. INTRODUCTION
Excited-state molecular dynamics (MD) simulations generally involve more than one potential energy surface (PES), since they need to describe both the relaxation and rearrangement processes in the initially populated state as well as the subsequent decay to the electronic ground state, possibly via other excited states.Of crucial importance for such internal conversions are conical intersections where different electronic states become degenerate (enabling radiationless transitions).To treat the dynamics around conical intersections properly, it is essential to apply an electronic structure method that offers a balanced description of all states involved.Configuration interaction (CI) methods and their multireference (MR) variants satisfy this condition, and they are thus in principle suitable for excited-state dynamics.
2][3][4] In this approach, the molecule is always in one particular state, but stochastic hops between states are possible (for example, according to the fewest switches criterion). 5,6The nuclear motion is propagated classically on the current PES, while the electronic properties are computed quantum-mechanically on-the-fly at every nuclear configuration.For the latter purpose, we use the semiempirical OM2/MRCI method [7][8][9] to determine the required electronic energies, gradients, and nonadiabatic couplings. 10This method is fast and has been applied a a) Electronic mail: spoerkel@kofo.mpg.deb) Electronic mail: thiel@kofo.mpg.denumber of medium-sized organic chromophores 11 (for some recent examples, see Refs.12-17).
In TSH simulations at the CI level, full CI calculations are normally not feasible for any but the very smallest systems, because the inclusion of all possible excitations will generally create a huge number of configurations that cannot be handled in practice.Therefore one usually defines an active space of orbitals and allows only certain types of excitations from pre-selected reference configurations.Only the configurations that can be created within these restrictions are included in the actual multi-reference CI expansion.For OM2/MRCI, we usually choose an active space of around 10 orbitals, with a few chemically relevant reference configurations, and include single and double (MRCISD) excitations.
In this article, we address technical problems that can arise in TSH simulations with OM2/MRCI and offer a simple solution, namely, the use of adaptive time steps.The key idea is to replace the current time step recursively by smaller steps until the problem disappears.The standard time step is recovered as soon as the problematic region is left.This correction procedure is local and easy to implement, and it increases the computational effort only when necessary.

A. Conventional surface hopping
During TSH simulations with active-space CI treatments, the orbitals may change their character and order as a result of the movement of the atoms.This can cause problems in the definition of the active space.For example, if the active space is comprised of the m highest occupied and n lowest unoccupied orbitals (m + n frontier orbitals selected by energy), its composition will change when the order of orbitals is switched such that a previously inactive orbital becomes active (i.e., one of the m + n frontier orbitals).In this case, there will be a sudden discontinuous change in the OM2/MRCISD energy and also in the gradient, which will corrupt the dynamics.It is thus essential to keep track of the orbitals along the trajectory, by computing their overlap between successive steps, and to retain the qualitative composition of the active space (by keeping the "tracked" orbitals even if they move out of the frontier range, while disregarding "intruder" orbitals).Evidently, this tracking procedure will become less effective and more error prone with increasing step size.
Similar remarks apply to another commonly employed definition of the active space, in terms of orbital character.For example, in TSH simulations of essentially planar molecules, it may be advantageous to use an active space comprised of all occupied and unoccupied π orbitals.This will generally be a "sparse" active space, with inactive σ orbitals lying energetically between the active π orbitals.In this case, a change in the order of σ and π orbitals during TSH dynamics is more likely than in the case of an active space composed of frontier orbitals only, simply because there are more opportunities for orbital switches.In this situation, tracking the orbitals along the TSH trajectory is obviously even more essential to ensure that the composition of the active space is retained during the simulation.
A related problem concerns the calculation of the gradient during TSH simulations using active-space CI methods.For example, if an active orbital and an inactive occupied orbital become close in energy, they will start mixing so that the character of the corresponding active orbital starts to change.This will affect the interactions with the correlating single and double excitations in the CISD treatment and lead to inconsistent energies and gradients, basically because the orbital and CI coefficients are optimized separately rather than simultaneously.In OM2/MRCISD TSH simulations, the occurrence of such unphysical energies and gradients can indeed be observed and traced back to strong orbital rotations of the kind described.One obvious remedy would seem to be an extension of the active space locally at and around such problematic points.We have experimented with this option, which however only alleviates the problems without solving them. 18Again, issues with detrimental orbital rotations will generally become more severe with increasing step size in the TSH simulations.
Unphysical local gradients caused by orbital rotations are usually much too high.The integration of the Newton equations (over large time steps) will then lead to high velocities, unrealistic geometries, and an unphysical increase of the kinetic energy.As a consequence, TSH simulations with such faulty gradients will tend to violate the conservation of total energy in the case of NVE ensembles.
In our previous work, we took several precautions against the problems outlined above and adopted the following standard procedure.The default time step for the propagation of the nuclei was chosen to be rather small (0.1 fs).Orbital tracking was always activated with the default requirement that active orbitals had to retain their character to 70% in two consecutive steps, i.e., the scalar product of the orbital coefficients had to exceed 0.7 in absolute value.At the post-processing stage, every finished trajectory was checked for sudden energy jumps; it was disregarded if the energy difference between any two steps was higher than a threshold (default 10 kcal/mol).Overall energy conservation was not checked in our default procedure.

B. Surface hopping with adaptive time steps
As described in Subsection II A, the TSH simulations may sometimes arrive at particularly problematic geometries where the mixing and switching between active and inactive orbitals causes abrupt changes in the computed electronic properties from one step to the next one.These problems are local in the sense that they normally appear only in certain small parts of configuration space.We propose a simple yet effective solution that is based on the key idea to reduce the integration time step in these areas (i.e., only when necessary).This is done by adopting the procedure specified in Figure 1.
This procedure ensures that the simulation always reaches time t k+1 irrespective of the number of substeps needed between t k and t k+1 .The simulation results are saved to output files at the uniformly distributed sequence of times t i for later analysis.Figure 2 illustrates how the chosen algorithm may overcome problems between steps t 1 and t 3 .
In our procedure, the adaptive time step can be reduced iteratively up to a user-defined maximum number of iterations (l max ).In the lth iteration, it equals step ∆t 0 , The acceptance criteria and the maximum number of iterations are chosen by the user.For our test examples, we adopted rather stringent criteria requiring a relative total energy change of less than 0.01% and all orbital overlaps to be greater than 0.97 in absolute value.For the two selected test cases (see below), the first criterion requires the energy change to be smaller than 0.01 and 0.03 kcal/mol, respectively, in consecutive steps; we note that our code also supports direct input of the acceptable energy change (recommended option).Furthermore, we chose l max = 10 so that the step can be reduced at most by a factor of 1024 (please note that SCF convergence and other thresholds may need to be tightened accordingly).If problems persist even after l max tries, the whole simulation is stopped and discarded.

A. GFP chromophore model
As a representative first test case, we selected a GFP chromophore model, 4-(2-hydroxybenzylidene)-1Himidazol-5(4H)-one (OHBI), which had previously been investigated in our group by OM2/MRCI nonadiabatics dynamics. 12The orbital active space consisted of 12 electrons in 12 π orbitals and was thus of "sparse" character.One may therefore expect to encounter in this example many of the problems discussed above.
We exactly followed the previous computational protocol 12 but ran more trajectories for better statistical analysis.We started 1000 trajectories both for the old procedure and the new adaptive procedure.Some of the trajectories terminated during the simulation for technical reasons, for example, failure to achieve SCF convergence or unsuccessful orbital tracking.Others were rejected at the post-processing stage, because of unphysically high energy jumps in two consecutive steps (more than 10 kcal/mol, see above).
In Table I, we present for both procedures the numbers of started, successful (correct), and unsuccessful trajectories.In the latter case, we distinguish between post-processing rejection, SCF convergence failure, and other numerical reasons (mainly orbital tracking).Obviously, the number of successful correct trajectories greatly increases from ca. 53% when using constant default time steps of 0.1 fs to ca. 95% with adaptive time steps.This is mainly because we ensure energy conservation between any two consecutive steps in the adaptive framework, which leads to much fewer rejections in the post-processing and also to a more stable simulation with less convergence and numerical problems (avoiding unphysically high gradients that generate unrealistic geometries in the next step).In our old procedure with constant time steps, the total energy may be much too high at the end of the simulation because errors between successive steps may accumulate; this may not be caught by our previous post-processing scripts that only check for an overly large energy difference between two steps (and not for overall energy conservation).
Figures 3 and 4 show the time evolution of the total energy for all trajectories and for the successful ones, with both procedures.Evidently, the total energy is not at all constant when using the constant default time step: even after the post-processing, the total energy increases significantly along most successful trajectories, on an average from 101.4 to 140.4 kcal/mol.When using the adaptive time step, only very few trajectories show a change in the total energy at all, and the average total energy remains essentially constant in the successful trajectories, with a small increase from 101.4 to 101.6 kcal/mol.In the old procedure, the observed increase in the total energy may cause artifacts in longer simulations, because the system may overcome barriers that are otherwise unsurmountable (leading to unrealistic photochemistry).
Figure 5 documents the actual number of steps taken along one particular (randomly chosen) trajectory in the adaptive scheme.For most of the time, the plot increases linearly  indicating that the default time step of 0.1 fs is adequate.
There are a few regions where the time step needs to be adapted because the acceptance criteria are not satisfied (see above).In these cases, the algorithm reduces the time step iteratively, additional substeps are taken, and the curve thus gets steeper.After leaving the problematic region, the time step adopts its default value again, and the curve becomes linear again.Obviously, in the chosen example, adaptive corrections are only needed on rare occasions, and the overall computation time increases only by about 10%.
Next we chose different time steps in the range of 0.05 fs up to 0.95 fs and ran simulations for 50 fs.The dependence of the number of adaptive corrections on the chosen time step is visualized in Figure 6.Shown are plots of the number of attempted steps (red and green arrows in Figure 2) and of accepted steps (green arrows in Figure 2) in the adaptive scheme; for comparison, we include a plot with the number of constant steps needed in the old procedure.
As expected, the adaptive algorithm does not need to intervene often for very small time steps (e.g., 0.05 fs) since the acceptance criteria are normally fulfilled.When the chosen time step becomes larger, the number of steps needed to reach 50 fs simulation time will generally tend to decrease.However, for larger time steps, the acceptance criteria will be fulfilled less often, so that the curves for the numbers of attempted and constant steps start to diverge.The plot for the number of attempted steps shows a minimum in the region of ∆t = 0.2 fs, where the time step is neither unnecessarily FIG. 5. Number of steps in the adaptive scheme for one particular trajectory.small nor large enough to require many adaptive corrections.When the chosen time step gets very large (e.g., beyond ∆t = 0.3 fs), the adaptive algorithm must reduce it very often so that the total number of attempted steps rises again.The data in Figure 6 indicate that the adaptive scheme is rather robust and that its computational demands do not depend much on the chosen step size.The adaptive procedure will always generate the steps that are required for getting stable results.The number of these steps seems to be minimal around ∆t = 0.2 fs but increases only slightly for larger default steps.In this sense, the adaptive scheme resembles a black-box method since the choice of time step does not have much impact.
Another distinction between the old procedure and the new adaptive procedure concerns the orbital tracking. 18,19The acceptance criterion for the mapping of the active orbitals has been tightened significantly: the absolute value of the scalar product of the orbital coefficients in two consecutive steps must exceed 0.70 in the old procedure compared with 0.97 in the new one.This is visualized for one particular (randomly chosen) trajectory in Figure 7, which contains these values for all accepted pairs of active orbitals.In the old procedure, there are quite a few smaller overlaps (in the range 0.7-0.9)which may lead to problematic behavior, whereas all overlaps are enforced to be close to 1 (above 0.97) in the new procedure.The use of a less stringent overlap criterion in the old procedure was a compromise between ensuring some degree of overlap and avoiding to have a huge percentage of failed trajectories.In the new procedure, a tight criterion can be imposed because the time step can be adapted sufficiently to satisfy the criterion also in problematic cases.
Finally, we address an obvious question: how much are our previous published results on the GFP chromophore 12 affected when using the new adaptive procedure?In our previous work, the decay from the S 1 to the S 0 state of the chromophore was studied by non-adiabatic dynamics simulations at the OM2/MRCI level, which yield the time evolution of the average state populations (and hence the lifetime of the excited state).The time evolutions are plotted FIG. 7. Overlaps in the mapping of the active orbitals during the tracking procedure for one particular trajectory of the GFP chromophore, for the old (top) and the new (bottom) procedure.
in Figure 8 as obtained from the old procedure and the new adaptive procedure.
It is reassuring that these plots are very similar and that the deactivation times are essentially identical (311 and 319 fs, FIG. 8. Time evolution of the average state occupations of the GFP chromophore obtained from all successful trajectories using constant and adaptive time steps.respectively, at the crossing points).While being technically superior, the new adaptive procedure thus gives essentially the same results as the old one in this particular test case.This may be related to the fact that the S 1 /S 0 deactivation of the GFP chromophore is an ultrafast process that is finished very early, before technical errors may accumulate.

B. Light-driven molecular motor
As second example, we consider a molecular motor studied recently in our group, 17 namely, a fluorinated N-alkylated indanylidene benzopyrrole (NAIBP), 3-[(2S)-2-fluoro-2-methyl-1-indanylidene]-1-methyl-2-methylindole (F-NAIBP), which was the technically most demanding system investigated by us so far, with the highest percentage of unsuccessful trajectories.The purpose of our previous work was the computational design of light-driven rotary molecular motors with improved quantum efficiency. 17OM2/MRCI TSH simulations gave quantum efficiencies of up to 47% for the four F-NAIBP conformers, roughly twice as high as those measured experimentally for the analogous allhydrocarbon motors. 17These results were based on more than 300 successful trajectories for each F-NAIBP conformer. 17n our previous report, we did not elaborate on the technical difficulties encountered in these OM2/MRCI TSH simulations.F-NAIBP has 41 atoms and is thus a rather large molecule, almost twice the size of the GFP model chromophore OHBI (22 atoms).The spectrum of valence orbital energies is therefore quite dense for F-NAIBP, significantly more so than in the case of OHBI, and the problems associated with the mixing and switching between active and inactive orbitals are more severe.In an attempt to cope with this problem, we chose a small time step of 0.05 fs (half the default value) in our previous work, 17 but still only a small fraction of trajectories finished successfully and fulfilled the continuity criteria.Hence we had to rerun most of the trajectories several times (with different random starting velocities) until they finished properly.TSH simulations with a constant time step are obviously very inefficient for F-NAIBP, which thus offers a stringent test for the proposed procedure with an adaptive time step.We chose the EP conformer 17 of F-NAIBP as the starting point for the present simulations.
For direct comparisons between the procedures with constant and adaptive time steps, we ran two sets of excitedstate dynamics with 600 trajectories each.In contrast to our previous work, we included all orbitals with π character in the active space of the MRCI treatment (8 electrons in 9 orbitals); previously we had left out one virtual orbital with π character that caused particularly severe problems in the TSH simulations. 17The default time step was presently chosen to be 0.1 fs.In the runs with constant time steps, trajectories were discarded if the overlap in the orbital mapping was below 0.70 or if there was an energy change of more than 10 kcal/mol between successive steps.When using adaptive time steps, the overlap in the orbital mapping was required to exceed 0.97 and the maximum relative change of total energy had to remain below 0.01% between successive steps.
The numbers of started and finished trajectories are listed in Table II.In the simulations with constant time steps, only ca.25% of the trajectories finish successfully and survive all checks.The number of successful trajectories is greatly increased to ca. 96% by the use of adaptive time steps; in this case, not a single trajectory is rejected in the post-processing checks on energy conservation.Figures 9 and 10 show the total energies along all trajectories for F-NAIBP.The plots are similar to those for OHBI (Figures 3 and 4), although after post-processing the results for constant and adaptive time steps differ less than those in the case of OHBI.The average total energy increases from 277.9 to 283.1 kcal/mol with constant time steps and  from 277.6 to 281.1 kcal/mol with adaptive time steps.Unlike OHBI, the total energy of F-NAIBP is thus not fully conserved even when using adaptive time steps, indicating the limitations of this procedure in difficult cases.
The average state occupations of F-NAIBP are shown in Figure 11.In contrast to the case of OHBI, there is a distinct difference in the computed lifetimes of the excited state: half of the trajectories have decayed to the ground state after 379 (452) fs for constant (adaptive) time steps.Apparently the use of adaptive time steps may allow trajectories to stay longer in the excited state by surviving visits to problematic regions, where trajectories are discarded in runs with constant time steps.
The central target property of our previous study 17 was the quantum efficiency of the photoinduced rotation of molecular motors, which is defined in terms of the dihedral angle of the central C= =C bond at the end of each trajectory (180 • at the outset).A final dihedral angle of less than 90 • is taken to indicate a successful rotation.The number of such "rotated" trajectories divided by the total number of trajectories gives the quantum efficiency.The present simulations yield values of 0.45 and 0.43 for trajectories with constant and adaptive FIG.11.Time evolution of the average state occupations of F-NAIBP obtained from all successful trajectories using constant and adaptive time steps.time steps, respectively.These are close to each other and to the value published previously for the EP conformer of F-NAIBP (0.47). 17 The technical improvements through the use of adaptive time steps thus lead to a notable increase in the lifetime of the excited state (by 19%) but have only a minor effect on the predicted quantum efficiency, thus corroborating the main qualitative conclusions from our previous study. 17

IV. DISCUSSION
In the proposed approach, adaptive time steps are introduced to overcome technical problems during the required on-the-fly electronic structure calculations, which are caused by the use of an active-space MRCI treatment that is, in principle, well suited for describing excited states and conical intersections.Of course, there are other quantum-chemical methods that do not suffer from such technical problems but face other limitations.One such example is linear-response density functional theory, which is widely used for excitedstate dynamics [1][2][3][4] but does not give the correct topology for conical intersections between the electronic ground state and the first singlet excited state. 20Because of their general applicability and solid theoretical foundation, MRCI methods are expected to remain a preferred choice in excited-state dynamics simulations, and technical advances such as the proposed adaptive time step protocol should therefore be of considerable practical value.
We are not aware of any other work in the TSH field that utilizes adaptive time steps for the purposes addressed in this article.Of course, adaptive protocols have proven to be valuable for other purposes and tasks in TSH simulations.For example, adaptive integration schemes have been developed for more accurate calculations of hopping probabilities, 21 and the fewest-switches algorithm has been modified by incorporating quantum uncertainty into the hopping times of classically forbidden hops. 22Recently, the occurrence of trivial crossings between noninteracting or weakly interacting states (e.g., in extended systems) has been recognized as a particular challenge for standard TSH simulations. 23,24In such cases, the nonadiabatic couplings peak sharply (being strongly localized near the exact crossing points) so that conventional fewest-switches TSH simulations require very small time steps and thus become very costly; in the noninteracting limit, the couplings become a delta function, and numerical simulations with a finite time step are no longer feasible. 25ence, these problems cannot be resolved by straightforward adaptation of the time steps, but need to be tackled by different strategies, including self-consistent fewest-switches surface hopping, 24 flexible surface hopping, 26 global flux surface hopping, 27 and local diabatization. 28,29While these novel schemes aim at improvements of the basic TSH formalism, our proposed adaptive time step protocol retains the standard TSH framework and aims at making TSH simulations with active-space MRCI treatments more robust and reliable.

V. CONCLUSIONS
The use of an adaptive time step can greatly improve the stability and accuracy of TSH simulations on potential surfaces computed with active-space CI methods.Problems arising from the mixing and switching between active and inactive orbitals can be controlled by sufficient reduction of the time step in such regions.The proposed algorithm supports energy conservation and successful orbital tracking also in problematic regions.It does not change the underlying methodology and increases the computational effort only slightly since smaller time steps are only applied locally when necessary.It is easily implemented and does not cause any changes in the post-processing analysis.The adaptive algorithm greatly lowers the number of trajectories that have to be rejected during post-processing.
Test calculations on a GFP chromophore model show that the adaptive scheme produces TSH trajectories with a much higher success rate (95%) than the previously used scheme (53%).Moreover, it ensures that the trajectories are technically sound in terms of energy conservation and orbital tracking.In the case of this GFP chromophore, these technical advances have no impact on the results of the TSH simulations: the previous protocol with constant time steps and the adaptive scheme yields essentially the same results for the time evolution of the average state occupations and for the deactivation times.
TSH simulations of the molecular motor F-NAIBP show an even bigger improvement in the success rate, which could be increased from 25% to 96%.In this case, trajectories that are kept alive through the use of adaptive time steps may apparently sample the excited-state PES for a longer time, which results in an average deactivation time that is about 70 fs larger than that obtained with constant time steps.On the other hand, the quantum efficiency of the photoinduced internal rotation of the molecular motor is hardly affected.
In summary, the use of adaptive time steps in TSH simulations with active-space CI methods is strongly recommended, because it helps to resolve technical problems, increases the efficiency and robustness of the simulations, and allows a more complete and realistic sampling of the PES.

FIG. 3 .
FIG.3.Time evolution of the total energy of the GFP chromophore with constant time step, for all trajectories (top) and for the successful trajectories (bottom).Shown in red: trajectories rejected in the post-processing.Green lines: average total energy and the associated standard deviation.

FIG. 4 .
FIG.4.Time evolution of the total energy of the GFP chromophore with adaptive time step, for all trajectories (top) and for the successful trajectories (bottom).Shown in red: trajectories rejected in the post-processing.Green lines: average total energy and the associated standard deviation.

194108- 5 LFIG. 6 .
FIG.6.Number of steps needed to reach 50 fs simulation time when using constant and adaptive steps.

FIG. 9 .
FIG. 9. Time evolution of the total energy of F-NAIBP with constant time step, for all trajectories (top) and for the successful trajectories (bottom).Shown in red: trajectories rejected in the post-processing.Green lines: average total energy and the associated standard deviation.

194108- 7 L
FIG.10.Time evolution of the total energy of F-NAIBP with adaptive time step, for all trajectories (top) and for the successful trajectories (bottom).Shown in red: trajectories rejected in the post-processing.Green lines: average total energy and the associated standard deviation.

TABLE I .
Number of different kinds of trajectories of the GFP chromophore for constant and adaptive time steps.

TABLE II .
Number of different kinds of trajectories for F-NAIBP when using constant and adaptive time steps.