# Fewest switches surface hopping with Baeck-An couplings

Mariana T. Do Casal, Josene M Toldo, Max Pinheiro Jr, Mario Barbatti

► To cite this version:

Mariana T. Do Casal, Josene M Toldo, Max Pinheiro Jr, Mario Barbatti. Fewest switches surface hopping with Baeck-An couplings. Open Research Europe, 2021, 1, 10.12688/openreseurope.13624.1. hal-03249607

HAL Id: hal-03249607

https://hal.science/hal-03249607

Submitted on 8 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

V1

First published: 17 May 2021, 1:49

https://doi.org/10.12688/openreseurope.13624.1

Latest published: 17 May 2021, 1:49

https://doi.org/10.12688/openreseurope.13624.1

## Abstract

In the Baeck-An (BA) approximation, first-order nonadiabatic coupling vectors are given in terms of adiabatic energy gaps and the second derivative of the gaps with respect to the coupling coordinate. In this paper, a time-dependent (TD) BA approximation is derived, where the couplings are computed from the energy gaps and their second time-derivatives. TD-BA couplings can be directly used in fewest switches surface hopping, enabling nonadiabatic dynamics with any electronic structure methods able to provide excitation energies and energy gradients. Test results of surface hopping with TD-BA couplings for ethylene and fulvene show that the TD-BA approximation delivers a qualitatively correct picture of the dynamics and a semiquantitative agreement with reference data computed with exact couplings. Nevertheless, TD-BA does not perform well in situations conjugating strong couplings and small velocities. Considered the uncertainties in the method, TD-BA couplings could be a competitive approach for inexpensive, exploratory dynamics with a small trajectories ensemble. We also assessed the potential use of TD-BA couplings for surface hopping dynamics with time-dependent density functional theory (TDDFT), but the results are not encouraging due to singlet instabilities near the crossing seam with the ground state.

## Keywords

Computational chemistry, Ultrafast dynamics, Nonadiabatic dynamics, Surface hopping, Fewest Switches, Ethylene, Fulvene

<figure><img src="image_3.png" /></figure>

This article is included in the Excellent Science gateway.

## Open Peer Review

### Reviewer Status

<figure><img src="image_4.png" /></figure>

### Invited Reviewers

1

2

3

version 1

17 May 2021

<figure><img src="image_5.png" /></figure>

report

<figure><img src="image_6.png" /></figure>

report

<figure><img src="image_7.png" /></figure>

report

1. **Sebastian Fernadez-Alberti**, Universidad Nacional de Quilmes, Bernal, Argentina

2. **Julia Westermayr** iD, University of Warwick, Coventry, UK

3. **Lluís Blancafort** iD, University of Girona, Girona, Spain

Any reports and responses or comments on the article can be found at the end of the article.
**Corresponding author:** Mario Barbatti (mario.barbatti@univ-amu.fr)

**Author roles:** **T. do Casal M:** Investigation, Validation, Writing - Original Draft Preparation; **Toldo JM:** Supervision, Validation, Writing - Original Draft Preparation;

**Pinheiro Jr M:** Data Curation, Methodology, Software, Validation, Writing - Original Draft Preparation;

**Barbatti M:** Conceptualization, Data Curation, Formal Analysis, Funding Acquisition, Investigation, Methodology, Project Administration, Software, Supervision, Visualization, Writing - Original Draft Preparation

**Competing interests:** No competing interests were disclosed.

**Grant information:** MP and MB were financially supported by the European Union's Horizon 2020 research and innovation programme under the grant agreement No 832237 (project SubNano). MTC, JMT, and MB were financially supported by the European Union's Horizon 2020 research and innovation programme under the grant agreement No 828753 (project BoostCrop). This work used computational resources of the project Equip@Meso (ANR-10-EQPX-29-01) funded by the French Government "Investissements d'Avenir" program.

**Copyright:** © 2021 T. do Casal M et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**How to cite this article:** T. do Casal M, Toldo JM, Pinheiro Jr M and Barbatti M. Fewest switches surface hopping with Baeck-An couplings [version 1; peer review: 3 approved] Open Research Europe 2021, 1:49 https://doi.org/10.12688/openreseurope.13624.1

**First published:** 17 May 2021, 1:49 https://doi.org/10.12688/openreseurope.13624.1
# Plain language summary

When a molecule absorbs light, its electrons are excited. The molecular geometry distorts until the excess of energy dissipates as heat or is reemitted as light. It is possible to simulate using a computer how a molecule reacts to light excitation. Still, such dynamics simulations involve complex and slow quantum-mechanical methods. In this work, we present a new approach to compute one of the quantities needed for the simulations, the nonadiabatic couplings. Based on a new theory proposed by Baeck and An, our method (we call it TD-BA) gets the couplings from electronic energies and their variations with time. Because these quantities are always available during dynamics simulations, TD-BA couplings are much simpler and faster to compute than usual couplings. We tested our TD-BA method for the simulation of two molecules, ethylene and fulvene. The results show that TD-BA does a good job predicting the dynamics of these molecules.

# I. Introduction

Surface hopping has become one of the main tools for nonadiabatic dynamics simulations of photoexcited molecules$^{1,2}$. Its popularity stems from its underlying independent trajectory and local approximations, requiring electronic quantities computed only at the nuclear geometry of classical trajectories. Moreover, surface hopping is usually done in adiabatic representation, avoiding cumbersome diabatization procedures usually required for wavepacket propagation. Thus, surface hopping is commonly implemented to work with on-the-fly calculations of electronic properties (energies, energy gradients, and couplings) obtained with many different quantum chemical methods$^{3-9}$, and it requires neither pre-calculated nor modeled global potential energy surfaces.

The primary constraint for interfacing surface hopping to a new quantum chemical method or program is to get the first-order nonadiabatic coupling vectors $\mathbf{h}_{JL} = \langle J | \frac{\partial}{\partial \mathbf{R}} L \rangle$ between adiabatic electronic states J and L. These quantities are not as widespread in quantum-chemical packages as excited-state energy gradients for instance, though the situation has improved in the last two decades$^{10-12}$. Alternatively, the coupling vectors can also be obtained from an approximation involving potential energy Hessian matrices$^{13,14}$. Moreover, there are surface hopping variants that do not even require explicit nonadiabatic coupling calculations, as is the case of the Zhu-Nakamura$^{15,16}$ and Belyaev-Lebedev$^{17-19}$ surface hopping approaches. The former estimates nonadiabatic events from minimum energy gaps and energy gradients around this minimum, while the latter does from minimum energy gaps and second time-derivative of the gap at the minimum gap. The popular fewest switches surface hopping (FSSH)$^{20}$ requires nonadiabatic coupling vectors. Nonetheless, because these vectors are always projected on the velocities within the formalism, they can be replaced by time-derivative couplings $\sigma_{JL} = \langle J | \frac{\partial}{\partial t} L \rangle = \mathbf{v} \cdot \mathbf{h}_{JL}$, which are much simpler to implement than nonadiabatic coupling vectors. In general, the calculation of time-derivative couplings follows the Hammes-Schiffer-Tully prescription$^{21}$, and they are obtained from wavefunction overlaps$^{22-24}$. Alternatively, the coupling

vector calculation can be avoided entirely in FSSH, by employing the local diabatization formalism$^{25}$, which is also based on wavefunction overlaps.

All these models possess algorithmic limitations for their universal adoption. The Hessian-based nonadiabatic coupling approximation is too costly for on-the-fly dynamics, although it has been used with machine-learning (ML) nonadiabatic dynamics$^{26}$. Zhu-Nakamura and Belyaev-Lebedev approaches require determination of the diabatic crossing point, which, translated into adiabatic representation, means that the hopping at time t depends not only on the previous timesteps but also on the future of the trajectory. (These methods are usually implemented in terms of a three-point approximation$^{15,17-19}$, evaluating the hopping at t from information computed at t - Δt, t, and t + Δt.) On the other hand, the wavefunction overlaps required by time-derivative coupling and local diabatization approaches boil down to nonorthogonal atomic-orbital overlap matrices, demanding lengthy algorithmic intervention whenever a new method or program needs to be interfaced to surface hopping.

Recently, Baeck and An$^{27}$ conjectured that nonadiabatic couplings near an energy crossing point could be estimated from the simple 1D model

$$h_{JL}(Q) \approx \frac{\operatorname{sgn}(\Delta E_{JL}(Q))}{2} \sqrt{\frac{1}{\Delta E_{JL}(Q)} \frac{\partial^2 \Delta E_{JL}(Q)}{\partial Q^2}} \quad Q \approx Q_c \quad (1)$$

involving only the energy gap $\Delta E_{JL} = E_J - E_L$ and the second derivative of this gap with respect to the coupling direction Q. Details of their model are surveyed in Section II. Baeck and An's numerical tests for different molecules have indeed shown an outstanding agreement between this simple approximation and the exact nonadiabatic coupling$^{27,28}$. Our goal in this paper is to expand and explore the Baeck-An (BA) model to be used in surface hopping dynamics.

Baeck and An pointed out two potential limitations of their model$^{27}$. First, it is strictly valid for molecular configurations coupling only two states. Second, the model is restricted to a single degree of freedom. When using their model for dynamics, there is nothing that we can do about the first limitation but to hope that two-state coupling dominates the nonadiabatic interactions. However, concerning the second limitation, it is fortunate that FSSH does not require the full nonadiabatic coupling vector but only its projection on the velocity direction. As we shall discuss, we can explore this feature to expand the 1D BA model to multidimensional dynamics, using a time-dependent (TD) BA model.

The TD-BA model is not supposed to replace the exact calculation of nonadiabatic couplings in surface hopping. It can be, however, an alternative approach to be used with electronic structure methods for which nonadiabatic couplings are still unavailable or are too computationally expensive to be computed on-the-fly. The TD-BA model can also be used as an inexpensive diagnostic of nonadiabatic interactions during exploratory dynamics, triggering higher-level calculations. Another
potential use of the TD-BA model is in ML. Several authors have reported difficulties in the ML-prediction of nonadiabatic couplings<sup>29–31</sup>. Their narrow functional shapes require too large training sets and active learning, being the main challenge for developing ML nonadiabatic dynamics. The TD-BA model, requiring only predicting energy gaps, may open a new avenue to overcome these problems.

In this work we show how TD-BA couplings can be seamlessly integrated into FSSH. We tested it by performing decoherence-corrected (DC) FSSH dynamics for ethylene and fulvene (Figure 1) using TD-BA couplings and comparing the results to DC-FSSH dynamics of these molecules using full nonadiabatic coupling vectors. Dynamics of both systems have been extensively studied before, and, recently, Ibele and Curchod<sup>32</sup> proposed that they may be considered molecular examples for two of the three basic 1D Tully models<sup>20</sup>. Ethylene is a crucial prototype<sup>32</sup> for internal conversion in various systems as diverse as protonated Schiff bases<sup>33</sup>, substituted polyenes<sup>34</sup>, and purine bases<sup>35</sup>. In turn, fulvene poses the challenge of S<sub>1</sub>↔S<sub>0</sub> recurrences at an extended crossing seam with both peaked and sloped conical intersections<sup>36</sup>. The results for both molecules are encouraging, showing that TD-BA-based dynamics delivers a semiquantitative agreement with reference datasets computed with exact couplings.

## II. The BA model for nonadiabatic couplings

Baeck and An have shown that the nonadiabatic coupling maximum for a two-states system can be approximated as<sup>27</sup>

$$h_{JL}(Q_c) \equiv \left\langle J \left| \frac{\partial L}{\partial Q} \right| Q_c \right\rangle = \frac{\operatorname{sgn}(\Delta E_{JL}(Q_c))}{2} \sqrt{\frac{1}{\Delta E_{JL}(Q_c)} \frac{\partial^2 \Delta E_{JL}}{\partial Q^2} \Big|_{Q_c}} \quad (2)$$

where Q is a one-dimensional coordinate coupling the adiabatic states J and L. Q_c is the geometry where the magnitude of the coupling is the largest. This expression arises from assuming that the nonadiabatic coupling can be represented as a Lorenz function of Q,

$$h_{JL}(Q) = \frac{\kappa \left( \frac{\Delta_c}{2} \right)}{2 \left( \left( \frac{\Delta_c}{2} \right)^2 + \kappa^2 (Q - Q_c)^2 \right)} \quad (3)$$

where κ, the coupling constant, and Δ<sub>c</sub>, the adiabatic energy gap at Q<sub>c</sub>, are constants. At Q<sub>c</sub>, the maximum nonadiabatic coupling is simply

<figure><img src="image_8.png" /><figcaption>Figure 1. Molecular structure of ethylene and fulvene.</figcaption></figure>

$$h_{JL}(Q_c) = \frac{\kappa}{\Delta_c} \quad (4)$$

Δ<sub>c</sub>, which by construction corresponds to twice the diabatic coupling, is given as the adiabatic energy gap between J and L at Q<sub>c</sub>.

$$\Delta_c \equiv \Delta E_{JL}(Q_c) \quad (5)$$

In turn, to get the coupling constant κ, Baeck and An explore the similarities between their model and the linear vibronic coupling (LVC) model<sup>37</sup>. While in LVC, Köppel, Gronki, and Mahapatra propose to compute κ from

$$\kappa = \left[ \frac{1}{8} \frac{\partial^2 \Delta E_{JL}^2}{\partial Q^2} \right]_{Q_c}^{1/2} \quad (6)$$

Baeck and An propose the expression

$$\begin{aligned} \kappa &= \left[ \frac{1}{8} \frac{\partial^2 \Delta E_{JL}^2}{\partial Q^2} - \frac{1}{4} \left( \frac{\partial \Delta E_{JL}}{\partial Q} \right)^2 \right]_{Q_c}^{1/2} \\ &= \frac{1}{2} \left[ \Delta E_{JL} \frac{\partial^2 \Delta E_{JL}}{\partial Q^2} \right]_{Q_c}^{1/2} \end{aligned} \quad (7)$$

which, compared to LVC, equivales to introduce an asymmetric linear term in Q in the expansion of ΔE<sub>JL</sub>. This term aims to make the model more flexible so that the minimum of ΔE<sub>JL</sub> does not necessarily coincide with Q<sub>c</sub>.

Finally, the nonadiabatic coupling maximum in Equation (2) is obtained by replacing κ given in Equation (7) into Equation (4).

Baeck and An also conjectured that the functional form of Equation (2) extended to an arbitrary point Q should describe the nonadiabatic coupling also in the region near Q<sub>c</sub>. This is the source of Equation (1).

## III. TD-BA couplings

To employ the BA model in surface hopping, we start by writing their expression as a vector in the unitary direction of Q

$$\mathbf{h}_{JL}(\mathbf{R}) \approx \frac{\operatorname{sgn}(\Delta E_{JL}(\mathbf{R}))}{2} \sqrt{\frac{1}{\Delta E_{JL}(\mathbf{R})} \frac{\partial^2 \Delta E_{JL}(\mathbf{R})}{\partial Q^2}} \hat{\mathbf{Q}} \quad \mathbf{R} \approx \mathbf{R}_c \quad (8)$$

where R is the full-dimensional molecular geometry, which has a crossing at R<sub>c</sub>. In the spirit of the Baeck-An 1D-approximation, Equation (8) assumes that the only non-null energy Hessian component is that one in the Q direction.

Using the chain rule, we can rewrite the Hessian as a time derivative

$$\frac{\partial^2 \Delta E_{JL}}{\partial Q^2} = \frac{1}{v_Q^2} \frac{d^2 \Delta E_{JL}}{dt^2} \quad (9)$$

where v<sub>Q</sub> is the projection of the nuclear velocity on Q̂. Replacing Equation (9) in Equation (8) yields the time-dependent version of the Baeck-An nonadiabatic coupling
$$ \mathbf{h}_{JL}(t) \approx \frac{\operatorname{sgn}(\Delta E_{JL}(t))}{2v_Q} \sqrt{\frac{1}{\Delta E_{JL}(t)}} \frac{d^2 \Delta E_{JL}(t)}{dt^2} \hat{\mathbf{Q}} \quad t \approx t_c \quad (10) $$

which implicitly assumes that $\hat{\mathbf{Q}}$ direction is constant. In the previous expression, $t_c$ is the time at which the molecule reaches $\mathbf{R}_c$.

In FSSH²⁰, the quantity of interest is the time-derivative coupling $\sigma_{JL} \equiv \langle J \left| \frac{d}{dt} L \right| J \rangle = \mathbf{v} \cdot \mathbf{h}_{JL}$ (see Section IV), which using Equation (10) simplifies to

$$ \sigma_{JL} \approx v_Q h_{JL} = \frac{\operatorname{sgn}(\Delta E_{JL})}{2} \sqrt{\frac{1}{\Delta E_{JL}}} \frac{d^2 \Delta E_{JL}}{dt^2} \quad (11) $$

Note that all references to $Q$ cancel out in Equation (11).

To enforce the condition $t \approx t_c$, we notice that whenever the molecule moves near a minimum of $\Delta E_{JL}$, the square-root argument in Equation (11) must be positive. Thus, we assume that the coupling is null for negative arguments, and the TD-BA nonadiabatic coupling becomes

$$ \sigma_{JL} \equiv \mathbf{v} \cdot \mathbf{h}_{JL} \approx \begin{cases} \frac{\operatorname{sgn}(\Delta E_{JL})}{2} \sqrt{\frac{1}{\Delta E_{JL}}} \frac{d^2 \Delta E_{JL}}{dt^2} & \text{if } \frac{1}{\Delta E_{JL}} \frac{d^2 \Delta E_{JL}}{dt^2} > 0 \\ 0 & \text{if } \frac{1}{\Delta E_{JL}} \frac{d^2 \Delta E_{JL}}{dt^2} \le 0 \end{cases} \quad (12) $$

Although the conditions in Equation (12) are not sufficient to detect $t \approx t_c$ (the square root argument may still be positive away from the minimum gap), they make a good job screening the relevant crossing regions during dynamics, as our tests show.

## IV. Using the TD-BA model in surface hopping

The FSSH hopping probability from $L$ to $J$ is computed as²⁰

$$ P_{L \to J}(t) = \max \left[ 0, \frac{-2\Delta\tau}{|c_L(t)|^2} \sigma_{LJ}(t) \operatorname{Re}(c_J(t)c_L^*(t)) \right] \quad (13) $$

where $\sigma_{LJ}$ is the time derivative coupling given by the TD-BA approximation in Equation (12), and the complex coefficients $c$ are the solution of the locally-approximated time-dependent electronic Schrödinger equation

$$ \frac{dc_J}{dt} = -\sum_K \left( \frac{i}{\hbar} E_J + \sigma_{JK} \right) c_K \quad (14) $$

$\Delta\tau$ in Equation (13) is the timestep for integrating Equation (14). It is, in general, much smaller than $\Delta t$, the timestep to integrate the classical equations of motion.

As usual, the FSSH over-coherence should be corrected and, here, we adopted the simplified decay of mixing³⁸.

Because the $\mathbf{h}$ direction is unknown, the velocity adjustment after hopping should be made in the direction of the linear momentum $\mathbf{p}$³⁹.

## V. Computational details and datasets

We have used ethylene and fulvene (Figure 1) dynamics to test the method. All calculations were done with the multiconfigu- rational self-consistent field (MCSCF) method state-averaged over three states for ethylene and two states for fulvene. Details on the computation level of ethylene are discussed in Ref. 39. In brief, the active space was composed of five separate complete active subspaces: the first subspace contained four electrons and four orbitals ($\pi$, $\sigma$, $\pi^*$, $\sigma^*$), and each of the other four subspaces contained two electrons in two orbitals ($\sigma$, $\sigma^*$). Thus, the active space can be denoted as [CAS(4,4) ⊕ 4×CAS(2,2)]. For fulvene, we adopted a stand- ard CAS(6,6) complete active space. The 6-31G(d) basis set⁴⁰ was used in both cases.

Dynamics was simulated with decoherence-corrected³⁸ fewest switches surface hopping²⁰ for a maximum of 200 fs for ethyl- ene and 60 fs for fulvene. The classical equation were integrated with $\Delta t = 0.1$ fs steps. The quantum equations were integrated with $\Delta \tau = 0.005$ fs steps, using interpolated electronic quan- tities between classical steps. Decoherence corrections were applied with the standard 0.1 au parameter³⁸.

For ethylene, one set of 500 trajectories was run with the TD-BA model. For the analysis, the TD-BA dataset was com- pared to surface hopping based on exact nonadiabatic coupling vectors with $\mathbf{h}$-adjusted and $\mathbf{p}$-adjusted velocities. These two datasets, also composed of 500 trajectories each, are the same ones discussed in Ref. 39. For fulvene, we run three data- sets of 200 trajectories each, one with the TD-BA model, and other two with exact nonadiabatic coupling vectors using $\mathbf{h}$-adjusted and $\mathbf{p}$-adjusted velocities.

For both molecules, the TD-BA model was applied with analytical second time derivatives from quadratic regres- sion ($\Delta T = 0.4$ fs) and $\delta\eta = 0.1$ au. The conditions asso- ciated to $\delta\epsilon$ and $\delta\bar{\omega}$ parameters were not applied. See Subsection VI.A for the definition and discussion of these four parameters.

The initial conditions were sampled from a harmonic oscilla- tor Wigner distribution of the nuclei. For ethylene, they were restricted to the 9.30 ± 0.25 eV excitation window³⁹. For ful- vene, the initial conditions were restricted to the 4.00 ± 0.34 eV window.

The MCSCF calculations were done with Columbus (ver- sion 7, 09-Oct-2020)⁴¹,⁴². Dynamics was done with Newton-X (version 2.2 build 15)⁷,⁴³ interfaced to Columbus. DC-FSSH based on TD-BA couplings with all options discussed in this paper is available in Newton-X version 2.2 build 15 or higher, to be used with any of the electronic structure methods and programs interfaced to Newton-X.

All datasets for ethylene and fulvene, with all Newton- X input and output files, are freely available (see Data availability)⁴⁴⁻⁴⁶.
Additional TD-BA DC-FSSH calculations were also done with linear-response time-dependent density functional theory (TDDFT) for fulvene and thiophene. Both sets of simulations were performed with the 6-31G(d,p) basis set using the CAM-B3LYP functional<sup>47</sup> for fulvene and ωB97XD functional<sup>48</sup> for thiophene. The classical equations of motion were integrated with $\Delta t = 0.1$ fs timestep for fulvene (maximum of 60 fs) and $\Delta t = 0.5$ fs for thiophene (maximum of 300 fs). The quantum equations were integrated with $\Delta \tau = \Delta t/20$, using interpolated electronic quantities between classical steps. The initial conditions were sampled from a harmonic oscillator Wigner distribution of the nuclei. For fulvene, they are the same used for complete active space self-consistent field (CASSCF). For thiophene, the initial conditions were restricted to the 6.2 ± 0.5 eV excitation window. The TD-BA model was applied with analytical second time derivatives from quadratic regression ($\Delta T = 0.4$ fs) and $\delta\eta = 0.1$ au. The conditions associated to $\delta\epsilon$ and $\epsilon\bar{\omega}$ parameters were not applied. Dynamics were done with Newton-X (version 2.2 build 15) interfaced to Turbomole 7.5<sup>49,50</sup>. GAMESS is a free and open source program that could be used in place of Turbomole to do the TDDFT calculations<sup>51</sup>.

The statistical analysis followed the protocol proposed in Ref. 39 to compare datasets. For ethylene, mean values and error bars (95% confidence interval) for seven observables (adiabatic-population and oscillator strength time constants and five structural channel yields) were computed from the trajectories. For fulvene, mean values and error bars were computed for six observables (first decay time constant, $S_1$ population at 20 fs, $S_1$ population at 60 fs, and three structural yields). These quantities were used to calculate the overlap score $\lambda_x$ for each observable and the mean overlap score $\Lambda^{(1,2)}$ for the average over all $\lambda_x^{39}$. $\lambda_x$ measures the probability that the observable predictions in two datasets overlap (for some confidence interval). In turn, $\Lambda^{(1,2)}$ is near one if the predictions of both datasets significantly overlap and zero if they do not.

# VI. Implementation and tests

## A. General aspects

The second time derivative in Equation (12) can be computed numerically. In our implementation, in the first two classical timesteps of the trajectory and after any hopping, the non-adiabatic coupling is assumed to be null. At the third step, the second time derivative is computed by finite differences with the backward $O(\Delta t)$ approximation

From the fourth timestep on, we use the backward $O(\Delta t^2)$ approximation

Alternatively, the second time derivative can be evaluated analytically after carrying out a quadratic regression over a

sequence of $N_p = \Delta T/\Delta t$ points accumulated over a period $\Delta T$ of a trajectory integrated with timestep $\Delta t$. After the regression, the predicted energy gap as a function of time is $\Delta E_{JL} \approx at^2 + bt + c$ and the second time derivative is

$$ \frac{d^2 \Delta E_{JL}}{dt^2} \approx 2a \qquad (17) $$

Again, the nonadiabatic coupling is assumed to be null in the first two timesteps of the trajectory and after a hopping, while data are accumulated to do the regression.

To reduce the instabilities in the derivatives caused by small energy discontinuities during dynamics (they are widespread in MCSCF propagation<sup>52</sup>), we adopted three “cleaning” conditions. Thus, in addition to the choice of second time derivative method (and the associated $\Delta T$ value in the case of using quadratic regression), these conditions introduce three parameters — $\delta\eta$, $\delta\epsilon$, and $\delta\bar{\omega}$ — whose values are discussed later in this section.

The first of such conditions checks the magnitude of the coupling variation in one timestep $\Delta t$. If the rate of this variation is larger than a parameter $\delta\eta$, the current coupling is discarded in favor of the coupling in the previous timestep:

$$ \bullet \text{ if } \frac{|\sigma_{JL}(t)| - |\sigma_{JL}(t - \Delta t)|}{\Delta t} > \delta\eta, \text{ then } \sigma_{JL}(t) = \sigma_{JL}(t - \Delta t) $$

The second condition eliminates large couplings at energy gaps that are too large. TD-BA couplings are computed only for energy gaps smaller than a threshold $\delta\epsilon$ and assumed to be zero otherwise:

$$ \bullet \text{ if } |\Delta E_{JL}(t)| > \Delta\epsilon, \text{ then } \sigma_{JL}(t) = 0 $$

The third condition avoids large couplings arising from points with large energy gaps but an anomalously large second time-derivative. The coupling is assumed to be null if the product $|\Delta E_{JL} \Delta E_{JL}''|$ is larger than an $\delta\bar{\omega}$ parameter:

$$ \bullet \text{ if } |\Delta E_{JL}(t) \Delta E_{JL}''(t)| > \delta\bar{\omega}, \text{ then } \sigma_{JL}(t) = 0 $$

Before doing any dynamics with the TD-BA couplings, we investigated the effect of using either method for evaluating second time derivatives and adopting the three cleaning conditions. To do so, we computed TD-BA couplings with different parameters setups for all points of the ethylene's h-adjusted dataset<sup>39</sup>. The results are summarized in Figure 2. This figure reports the root mean square deviation (RMSD) between the TD-BA coupling $\sigma_{JL}$ and the exact $\hat{\sigma}_{JL}$ computed for the h-adjusted dataset

$$ RMSD = \sqrt{\frac{1}{N} \sum_{k=1}^{N_{\text{Traj}}} \sum_{i} \left( |\sigma_{JL}^{k}(t_i)| - |\hat{\sigma}_{JL}^{k}(t_i)| \right)^2} \qquad (18) $$

where N is the total number of timesteps in all trajectories.

The RMSD in Figure 2 shows that couplings computed from analytical second derivatives from a quadratic regression are always more accurate than couplings computed from numerical second derivatives. The difference between the two
<figure><img src="image_9.png" /><figcaption>**Figure 2.** Root mean square deviation (RMSD) between time-dependent Baeck-An (TD-BA) and exact nonadiabatic couplings computed for the ethylene's *h*-adjusted dataset. (a) RMSD is given for different values of δη with analytical (from quadratic regression over Δ*T* = 0.4 fs) and numerical (via finite differences) second time derivatives. The energy gap threshold (δε) is 12 eV and δω = 0.01 au. (b) RMSD is given for different quadratic regression periods (Δ*T* for δη = 0.1 au, δε = 12 eV, and δω = 0.01 au. (c) RMSD is given for different energy gap threshold (δε) for δη = 0.1 au, Δ*T* = 0.4 fs, and δω = 0.01 au. (d) RMSD is given for different δω values for δη = 0.1 au, Δ*T* = 0.4 fs, and δε = 12 eV. Note the log scales in (a) and (d).</figcaption></figure>

procedures is particularly remarkable for δη larger than 1 au (Figure 2a). Nonetheless, the period Δ*T* over which the quadratic regression is done should be very short, smaller than 0.5 fs (Figure 2b). This means that for dynamics with large timesteps, using the quadratic regression may not be advantageous. The parameters δε (Figure 2c) and δω (Figure 2d) play only a minor role in the BA couplings' accuracy.

Density maps comparing the magnitudes of BA couplings to exact nonadiabatic couplings for all points of the *h*-adjusted dataset for ethylene are shown in Figure 3. As expected, most couplings are near zero. The largest couplings for each pair of states tend to fall on the diagonal, which implies a good agreement between TD-BA and the exact couplings. The order of magnitude is generally well represented, with the S₀-S₂ coupling magnitudes ten times smaller than those of the other two transitions. The TD-BA model tends to overestimate small-magnitude couplings, especially for the S₀-S₂ interactions. The model also shows some underestimation for large couplings in S₀-S₁ and S₁-S₂ interactions.

Figure 4 shows results comparing TD-BA and exact couplings for a single trajectory of the *h*-adjusted dataset. In this case, after starting in S₁, ethylene hops to S₂ after 12.5 fs and stays there until it comes back to S₁ at 48.1 fs. It finally returns to S₀ at 93 fs (Figure 4a). The TD-BA model delivers a semiquantitative agreement with the exact results for the coupling magnitudes between all three pairs of states (Figure 4b, c, and d). The model correctly predicts the time where each peak occurs, as well as gives the right shape (intensity and width) of the coupling peaks. However, the TD-BA model shows spurious spikes that even the cleaning conditions discussed above could not get rid of.

## B. Potential use of TD-BA in TDDFT dynamics

Linear-response TDDFT is the workhorse of excited-state calculations of large molecular systems⁵³. Nevertheless, there are some critical obstacles to simulate nonadiabatic dynamics based on this method, especially when dealing with the internal conversion to the ground state. Because of its underlying approximations⁵⁴, TDDFT not only fails in describing the
<figure><img src="image_10.png" /><figcaption>**Figure 3.** Density maps comparing the magnitudes of time-dependent Baeck-An (TD-BA) and exact nonadiabatic couplings for all points of each pair of states in the ethylene's *h*-adjusted dataset. TD-BA; smoothed over $\Delta T = 0.4$ fs, and with $\delta\eta = 1$ au, $\delta\epsilon = 12$ eV, and $\delta\omega = 0.01$ au. A low-level contour is shown in each graph to guide the eye. The maps were computed with 672,812 data points.</figcaption></figure>

state crossing between the ground and the first excited state, but also predicts the wrong dimensionality of the branching space around the crossing⁵⁵. Such handicaps have not precluded cautious use of TDDFT for studying internal conversion to the ground state. Indeed, thanks to Send and Furche¹⁰, $S_1/S_0$ nonadiabatic coupling vectors for TDDFT are now a common feature of several quantum chemistry software programs. Our own research protocol has been to run TDDFT dynamics until the energy $S_1/S_0$ gap drops below some threshold (about 0.1 eV), stop the trajectory propagation there, and assume that point to be the hopping point³⁴. Here, we were interested in testing whether the TD-BA model would allow going beyond this protocol, with actual hoppings to the ground state. To do so, we performed exploratory simulations with fulvene and thiophene.

The first challenge is to obtain an acceptable (TD) DFT description in regions close enough to the crossing region between the ground and excited state so that the hops between them can be computed. Therefore, we increased the size of the grid in the DFT calculations to the maximum value available in Turbomole 7.5 (grid size 7) as an attempt to avoid self-consistent field (SCF) convergence problems and negative excitation energies. Although this strategy successfully solved SCF convergence problems, our trajectories for fulvene were systematically terminated during the linear-response calculation of excitation energies due to singlet instabilities at small $S_0/S_1$ energy gaps. To further avoid singlet instabilities, we used the Tamm-Dancoff approximation (TDA)⁵⁶, as proposed by Tapavicza et al.⁵⁷ We computed 25 trajectories starting in the $S_1$ state, of which 22 did not show any hopping between $S_1$

<figure><img src="image_11.png" /><figcaption>**Figure 4.** Benchmark tests on one trajectory (TRAJ1) from the ethylene's *h*-adjusted dataset. (a) Potential energies as a function of time. Exact and time-dependent Baeck-An (TD-BA) coupling magnitudes during this trajectory are shown in (b) for $S_0-S_1$, (c) for $S_0-S_2$, and (d) for $S_1-S_2$.</figcaption></figure>
and S₀. In the remaining three trajectories, the hopping occurred at energy gaps of 0.41, 2.22, and 1.65 eV. This set of results suggests that the coupling between S₁ and S₀ is not well described for fulvene under the TDA approximation. Moreover, the S₁/S₀ energy gaps as a function of time were systematically too large (compared to those at CASSCF), indicating that the TDA description of the branching space was inadequate. Therefore, neither TDDFT nor TDA-DFT were suitable to be used for fulvene dynamics.

To test how system-dependent this failure was, we turned to thiophene. Surface hopping dynamics of thiophene with TDA was previously reported by Fazzi et al.,⁵⁸ who showed that TDA was able to describe the nonradiative relaxation process bringing thiophene to the S₁/S₀ crossing region. We computed 12 trajectories starting from the S₂ state. From those, six hopped to the ground state with reasonably small energy gaps. Nevertheless, five prematurely ended due to singlet instability at the crossing region. We tried to avoid the singlet instability by decreasing the time step to 0.1 fs, hoping to have a hopping before reaching the instability point. However, the instability persisted.

Therefore, although the TD-BA model delivers an alternative way to compute the nonadiabatic couplings between the ground and excited states, its application to nonadiabatic dynamics at TDDFT level is discouraging since it still depends on the good description of the relevant states close to the crossing region. We expect similar negative results for dynamics based on algebraic diagrammatic construction to second-order (ADC(2))⁵⁹, which tends to yield negative excitations near the S₁/S₀ crossing. Note, however, that TD-BA can still be helpful for both TDDFT and ADC(2) dynamics involving only internal conversion between excited states.

# VII. Dynamics results

## A. Ethylene

Ethylene dynamics has been extensively investigated with many theoretical approaches, including wavepacket dynamics⁶⁰, on-the-fly quasi-diabatic dynamics⁶¹, multiple spawning⁶²⁻⁶⁵ multiconfigurational Ehrenfest⁶⁶, and surface hopping³²,⁶⁷⁻⁷². The ultrafast nonadiabatic dynamics of ethylene is driven by torsional and pyramidalization modes until the molecule finds the S₁/S₀ crossing seam.

The adiabatic-state population evolution based on TD-BA couplings is shown in Figure 5 (top). The S₁ state is quickly depopulated toward S₀ within 100 fs. In the first 20 fs, S₁ also receives some of the population, but it quickly returns to the lower states. The S₁ lifetime fitted with an exponential decay is 107 ± 9 fs, where the margin of error is computed for a 95% confidence interval.

Population evolutions with TD-BA and exact couplings are also compared in Figure 5 (top). The results are encouraging, delivering a semiquantitative agreement between the two datasets. The S₁ lifetime in the *h*-adjusted dataset is 114 ± 10 fs, meaning that the results agree within the error bars. This level of variation is expected even within different sets of FSSH

<figure><img src="image_12.png" /><figcaption>**Figure 5.** Time evolution of time-dependent Baeck-An (TD-BA) decoherence-corrected fewest switches surface hopping (DC-FSSH) dynamics of ethylene. (Top) mean adiabatic populations of S₀, S₁, and S₂ (bottom) mean oscillator strength between the current and the ground states. "Exact" corresponds to DC-FSSH dynamics with exact couplings from the *h*-adjusted dataset.</figcaption></figure>

with exact couplings. For example, the lifetime of FSSH dynamics with exact couplings but adjusting velocities in the *p* direction is 120 ± 11 fs. Nevertheless, there are quantitative differences between TD-BA and the *h*-adjusted dataset. Although the TD-BA dynamics perfectly describes the S₂ population evolution, the S₁ and S₀ evolutions bear some significant differences between the two datasets, with the S₁ population after 50 fs always lower in the *h*-adjusted dataset.

The current state's mean oscillator strength is shown in Figure 5 (bottom) for TD-BA and the *h*-adjusted dataset. The time constant obtained from an exponential fitting is 14 ± 2 fs, in excellent agreement with the reference value, 15 ± 2 fs. TD-BA correctly describes the damped oscillations in ethylene dynamics, which have been previously predicted by wavepacket⁶⁰ and surface hopping⁶⁷,⁶⁸ simulations, and experimentally measured⁷³,⁷⁴.

One of the most significant differences between the TD-BA and the *h*-adjusted datasets is the fraction of molecules still in the excited states after 200 fs (Table 1). While the dynamics with exact couplings put this quantity at 9%, TD-BA dynamics has almost twice this amount, 16%. This disagreement is out of the margin of error of ±3%. Beyond this divergence, the ground-state-product' yields show much better agreement. At the end of the simulations, dynamics with TD-BA couplings predicts 38% of ground-state CH₂CH₂ structures,
<table><caption>**Table 1.** Mean value and error bars (95% confidence interval) of time constants, populations, and structural yields, computed for ethylene and fulvene dynamics. Ethylene's results for the h-adjusted (reference) and p-adjusted datasets of Ref. 39 are shown as well. $\lambda_x^h$ and $\lambda_x^p$ are the overlap scores for each observable computed between time-dependent Baeck-An (TD-BA) and the respective dataset.</caption><thead><tr><th>Observable</th><th>TD-BA</th><th>h dataset</th><th>λ<sub>x</sub><sup>h</sup></th><th>p dataset</th><th>λ<sub>p</sub><sup>h</sup></th></tr></thead><tbody><tr><td><strong>Ethylene</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Population time τ<sub>p</sub> (fs)</td><td>107±9</td><td>114±10</td><td>0.50</td><td>120±11</td><td>0.07</td></tr><tr><td>Osc. strength time τ<sub>o</sub> (fs)</td><td>14±2</td><td>15±2</td><td>0.73</td><td>17±2</td><td>0.02</td></tr><tr><td>S<sub>2</sub>, S<sub>1</sub> CH<sub>2</sub>CH<sub>2</sub> yield (%)</td><td>16±3</td><td>9±3</td><td>0.00</td><td>12±3</td><td>0.05</td></tr><tr><td>S<sub>0</sub> CH<sub>2</sub>CH<sub>2</sub> yield (%)</td><td>39±4</td><td>38±4</td><td>0.84</td><td>41±4</td><td>0.83</td></tr><tr><td>S<sub>0</sub> CH<sub>3</sub>CH yield (%)</td><td>14±3</td><td>18±3</td><td>0.13</td><td>15±3</td><td>0.94</td></tr><tr><td>S<sub>0</sub> C<sub>2</sub>H<sub>3</sub> + H yield (%)</td><td>28±4</td><td>31±4</td><td>0.60</td><td>27±4</td><td>0.83</td></tr><tr><td>S<sub>0</sub> C<sub>2</sub>H<sub>2</sub> + 2H yield (%)</td><td>1±1</td><td>3±2</td><td>0.22</td><td>3±1</td><td>0.01</td></tr><tr><td><strong>Fulvene</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>1<sup>st</sup> decay time τ<sub>t</sub> (fs)</td><td>11±1</td><td>9±1</td><td>0.00</td><td>9±1</td><td>0.00</td></tr><tr><td>S<sub>1</sub> population at 20 fs (%)</td><td>46±7</td><td>34±6</td><td>0.00</td><td>42±7</td><td>0.66</td></tr><tr><td>S<sub>1</sub> population at 60 fs (%)</td><td>12±5</td><td>3±3</td><td>0.00</td><td>6±3</td><td>0.02</td></tr><tr><td>Planar S<sub>1</sub>→S<sub>0</sub> hop (%)</td><td>85±5</td><td>87±5</td><td>0.81</td><td>84±5</td><td>0.93</td></tr><tr><td>Tw-Stretched S<sub>1</sub>→S<sub>0</sub> hop (%)</td><td>5±3</td><td>6±3</td><td>0.86</td><td>8±4</td><td>0.38</td></tr><tr><td>Tw-Shrunk S<sub>1</sub>→S<sub>0</sub> hop (%)</td><td>9±4</td><td>7±4</td><td>0.73</td><td>8±4</td><td>0.91</td></tr></tbody></table>

while dynamics with exact couplings, 39%. According to the TD-BA and reference sets, the expected amounts of ethyldene (CH<sub>3</sub>CH) are 14% and 18%, respectively. Single-H dissociation yields are 28% and 31%. Double-H dissociation yields are 1% and 3%. All these figures agree within their margin of errors.

The hopping counting and the mean energy gap at the hopping time for all transitions are shown in Table 2 for the three datasets. Concerning the mean gap, TD-BA is in excellent agreement with the other two sets for all transitions excepting S<sub>0</sub>→S<sub>2</sub>. This discrepancy is not statistically relevant given that there are only a few of such rare transitions in each set. The number of hoppings between S<sub>1</sub> and S<sub>0</sub> (both directions) also agrees well between TD-BA and the other two sets. Nevertheless, the number of hoppings between S<sub>1</sub> and S<sub>2</sub> (both directions) is twice as large in TD-BA as in the other sets. With exact couplings, 43% of trajectories populated S<sub>2</sub>. These trajectories returned to S<sub>1</sub> afterward and, in general, remained there until they hopped to S<sub>0</sub>. With TD-BA couplings, the fraction of trajectories populating S<sub>2</sub> was higher, 63%. These trajectories usually hopped forth and back between S<sub>1</sub> and S<sub>2</sub> twice before converting to S<sub>0</sub>.

We can understand the excess of S<sub>1</sub>-S<sub>2</sub> hoppings in the TD-BA approach by analyzing a single trajectory computed with

exact and TD-BA couplings, as illustrated in Figure 6. The comparison is made to a p-adjusted dataset trajectory, limiting the differences to the coupling itself, without the influence of the velocity adjustment after hopping.

In the initial stage of the S<sub>1</sub> dynamics, the CC stretching strongly stabilizes the zwitterionic state S<sub>2</sub> (π*(2)). As shown in Figure 6a, the S<sub>1</sub> and S<sub>2</sub> states approach each other every approximately 15 fs. S<sub>2</sub> destabilizes again when the CC bond shrinks due to the oscillatory motion. This initial oscillation of the S<sub>1</sub>-S<sub>2</sub> gap is a common feature of all trajectories. When dynamics is propagated with TD-BA couplings, the hopping probability from S<sub>1</sub> (ππ*) to S<sub>2</sub> (π*(2)) at the stretched structure is larger than when propagating with exact couplings. Thus, the TD-BA trajectories hop to S<sub>2</sub> about twice more often than the reference trajectories.

The reason for the larger probability is connected to the surface topography. With the exact coupling, the projection h on v is modulated by the velocity magnitude. When the CC stretching is maximum, the velocity in this internal coordinate is null, and v.h is null as well, causing the zero at 15 fs that we can see in the orange curve (p-adjusted dataset) in Figure 6c. TD-BA coupling fails to capture this dependence on |v|, yielding a larger coupling peak at the same time. Nevertheless, this discrepancy in the probability does not play a
<table><caption>**Table 2.** Number of hoppings (count) and hopping energy-gap absolute mean value and standard deviation for the three sets of trajectories of ethylene and fulvene.</caption><thead><tr><th rowspan="2"></th><th colspan="3">TD-BA set</th><th colspan="3">h set</th><th colspan="3">p set</th></tr><tr><th>Count</th><th>Mean (eV)</th><th>St. Dev. (eV)</th><th>Count</th><th>Mean (eV)</th><th>St. Dev. (eV)</th><th>Count</th><th>Mean (eV)</th><th>St. Dev. (eV)</th></tr></thead><tbody><tr><td><strong>Ethylene</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>S<sub>2</sub>→S<sub>1</sub></td><td>520</td><td>0.35</td><td>0.22</td><td>259</td><td>0.39</td><td>0.28</td><td>252</td><td>0.43</td><td>0.32</td></tr><tr><td>S<sub>2</sub>→S<sub>0</sub></td><td>38</td><td>4.63</td><td>0.42</td><td>18</td><td>4.93</td><td>0.47</td><td>20</td><td>4.91</td><td>0.97</td></tr><tr><td>S<sub>1</sub>→S<sub>2</sub></td><td>563</td><td>0.61</td><td>0.45</td><td>281</td><td>0.63</td><td>0.43</td><td>285</td><td>0.62</td><td>0.42</td></tr><tr><td>S<sub>1</sub>→S<sub>0</sub></td><td>546</td><td>0.51</td><td>0.63</td><td>618</td><td>0.65</td><td>0.88</td><td>613</td><td>0.61</td><td>0.85</td></tr><tr><td>S<sub>0</sub>→S<sub>2</sub></td><td>4</td><td>3.91</td><td>0.54</td><td>2</td><td>2.12</td><td>0.29</td><td>6</td><td>1.88</td><td>0.19</td></tr><tr><td>S<sub>0</sub>→S<sub>1</sub></td><td>160</td><td>0.71</td><td>0.81</td><td>181</td><td>0.76</td><td>0.86</td><td>187</td><td>0.70</td><td>0.82</td></tr><tr><td><strong>Fulvene</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>S<sub>1</sub>→S<sub>0</sub></td><td>273</td><td>0.25</td><td>0.27</td><td>273</td><td>0.28</td><td>0.28</td><td>278</td><td>0.25</td><td>0.25</td></tr><tr><td>S<sub>0</sub>→S<sub>1</sub></td><td>92</td><td>0.57</td><td>0.56</td><td>82</td><td>0.56</td><td>0.40</td><td>91</td><td>0.47</td><td>0.32</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 6.** Characterization of the first 50 fs of a single trajectory of ethylene. Potential energies of the S<sub>1</sub> and S<sub>2</sub> states as a function of time with (a) exact p-adjusted dataset and (b) time-dependent Baeck-An (TD-BA) couplings. The dots indicate the current state every 0.5 fs. Violet dots indicate ππ* diabatic character, and green dots indicate π*(2) character. The absolute value of the exact and TD-BA v.h projection computed for the reference trajectory is shown in (c). The potential energy profile of S<sub>1</sub> and S<sub>2</sub> along the CC distance, keeping the other coordinates at the S<sub>0</sub> minimum, is given in (d).</figcaption></figure>
significant role in the dynamics — at least not for ethylene. After populating S₂, the molecule quickly returns to S₁ within 5 fs or less, and the population evolution is not affected, as we have seen in Figure 5.

Table 1 gives the overlap score for all time constants and structural yields discussed above. When the TD-BA dataset is compared to the h-adjusted dataset, the scores $\lambda_x^h$ are excellent for the S₀ CH₂CH₂ yield (0.8) and oscillator strength time constant (0.7). For the population time constant (0.5) and single-H dissociation (0.6), the scores are fine. They are poor for the less populated channels, ethylidene (0.1) and double-H dissociation (0.2). The results do not agree at all for excited-state CH₂CH₂ (0.0).

When the TD-BA dataset is compared to the p-adjusted dataset ($\lambda_x^p$), the lifetime agreement deteriorates (0.07 score for population and 0.02 for oscillator strength), but most of the structural figures improve (Table 1). There is an almost perfect agreement for ground-state CH₂CH₂, CH₃CH₂, and single H dissociation yields between the two sets, with overlap scores above 0.8. This means that part of the divergence to the h-adjusted dataset is caused not by the TD-BA couplings but by the p-adjustment of the velocities, which directly impacts the fragmentation process.

The statistical analysis in this paper follows the methods and protocol defined proposed in Ref. 39 Thus, the current results can be directly compared to those previously published datasets. Table 3 compiles a list of mean overlap scores $\Lambda^{(1,2)}$, comparing different FSSH setups to the h-adjusted reference dataset. The idea is that the implementation of new methods and testing of different features can always be analyzed in the same way to build a hierarchical comparison of procedures. We see that TD-BA yields $\Lambda^{(1,2)}$ of 0.43 at the same level as that of the p-adjusted dataset (0.41). This is equivalent to say that the results of ethylene DC-FSSH dynamics with TD-BA couplings and with exact couplings adjusting the velocity in the p direction are statistically equivalent for 500 trajectories

and 95% confidence interval. Nevertheless, this statement does not mean that TD-BA results should be accepted without questioning. The average overlap score varies between 0 (complete disagreement with reference data) and 1 (perfect agreement). An average overlap score of 0.43 means that significant differences may be expected in comparison to the reference data. However, these differences are not worse than when we use exact couplings and do velocity adjustment in the p direction. As we discussed in the previous paper³⁹, if we reduce the number of trajectories in the ensemble, the margins of error grow, and all overlap scores improve. For a small but typical ensemble of 100 trajectories, the mean overlap score for TD-BA dynamics is about 0.7, which is an acceptable agreement with the reference data. Therefore, TD-BA dynamics seems to be a promising approach for exploratory dynamics with few trajectories.

## B. Fulvene

Fulvene nonadiabatic dynamics has been discussed before based on different simulation approaches, including wavepacket propagation in reduced-dimensionality⁷⁵⁻⁷⁹, full-dimensionality direct-surface wavepacket propagation⁸⁰,⁸¹, multiple spawning³², and surface hopping³²,⁸². Fulvene has an extended S₁/S₀ crossing seam from planar to 90° twisted structures³⁶,⁸²,⁸³. The central feature of its photophysics is the ultrafast sub-100 fs decay to the ground state, with periodic returns to S₁ and internal conversion spread over the entire crossing seam.

For the analysis that follows, it is helpful to define the mean torsional angle as the average of the absolute values of the four dihedral angles around the C-CH₂ bond, which in degrees is:

$$ \phi_{\text{C-CH}_2} = \frac{1}{4} \left[ |\phi_{\text{CC-CH}}^{\text{cis}_1}| + |\phi_{\text{CC-CH}}^{\text{cis}_2}| + |180^{\circ} - \phi_{\text{CC-CH}}^{\text{trans}_1}| + |180^{\circ} - \phi_{\text{CC-CH}}^{\text{trans}_2}| \right] \quad (19) $$

To compare the results of different datasets, we have defined three regions in the C-CH₂ torsion and distance space. The first region (Planar) includes near planar structures with a mean torsional angle smaller than 30°. The second region (Twisted-stretched) includes twisted structures (≥ 30°) with a C-CH₂ distance larger than 1.55 Å. The third region (Twisted-shrunk) includes twisted structures with C-CH₂ distances smaller than 1.55 Å. The count of S₁→S₀ hoppings in each region is given in Table 1. These three regions are indicated in Figure 7a-c.

The three datasets show that same general behavior. After the excitation, fulvene returns to the ground state within 10 fs, with S₁ repopulations (recurrences) at 20, 40, and 60 fs (Figure 7d). Internal conversion to S₀ can occur at different points of the crossing seam. As shown in Figure 7 (a-c), the first wave of S₁→S₀ hoppings occurs in the Planar region, in the first 20 fs after excitation. In the h-adjusted dataset (Figure 7b), the mean value of this angle is 13° (±7°), and the mean C-CH₂ distance is 1.58 Å (±0.06 Å). For comparison, the crossing seam has a maximum at 0° and 1.58 Å, implying that the first wave of S₁→S₀ hoppings occurs near this

<table><caption>**Table 3.** Mean overlap score $\Lambda^{(1,2)}$ computed for different fewest switches surface hopping (FSSH) setups. Ethylene: 500 trajectories in each dataset and 95% confidence interval averaged over seven observables. Fulvene: 200 trajectories in each dataset and 95% confidence interval averaged over six observables. The h- and p-adjusted datasets of ethylene are from Ref. 39.</caption><thead><tr><th rowspan="2">Dataset</th><th colspan="2">Λ<sup>(1,2)</sup></th></tr><tr><th>Ethylene</th><th>Fulvene</th></tr></thead><tbody><tr><td>h-adjusted (reference)</td><td>1</td><td>1</td></tr><tr><td>p-adjusted</td><td>0.41</td><td>0.59</td></tr><tr><td>TD-BA</td><td>0.43</td><td>0.40</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure 7. Characterization of fulvene dynamics.** C-CH₂ mean dihedral angle and the CC distance at the S₁→S₀ hopping for the (a) time-dependent Baeck-An (TD-BA), (b) *h*-adjusted, and (c) *p*-adjusted datasets. The colors indicate the hopping time. The dotted circle shows the Franck-Condon (FC) position, and the crosses show the positions of the planar conical intersection (0°), the minimum on the crossing seam (63°), and the twisted conical intersection (90°). The dashed curve indicates the minimum on the crossing seam optimized with constrained angular values. The vertical and horizontal lines split the space into three regions, Planar (φ<sub>C-CH₂</sub> < 30°), Twisted-stretched (φ<sub>C-CH₂</sub> ≥ 30°, d<sub>C-CH₂</sub> ≥ 1.55 Å), and Twisted-shrunk (φ<sub>C-CH₂</sub> ≥ 30°, d<sub>C-CH₂</sub> < 1.55 Å). Panel (d) shows the S₁ population as a function of time for the three datasets.</figcaption></figure>

maximum. The second S₁→S₀ hopping wave happens between 20 and 40 fs at stretched C-CH₂ distances (> 1.55 Å). In this second wave, the structure may be either planar or twisted. A third S₁→S₀ hopping wave happens between 40 and 60 fs, mostly at the Twisted-shrunk region. This third wave occurs near the minimum of the crossing seam, located at 63° and 1.48 Å.

To proceed with the statistical analysis, including the calculation of overlap scores, the S₁ population error bar (95% confidence interval) was estimated with bootstrap through 10,000 repetitions of the 200 trajectories of each dataset.

The overlap scores are above 0.7 for the number of S₁→S₀ hoppings in each of the three regions of the torsion × distance space (Table 1). It means that the TD-BA model quantitatively agrees with both **h**- and **p**-adjusted datasets for predicting these observables. The only exception is the divergence in the number of Twisted-stretched hoppings between TD-BA and the **p**-adjusted dataset (0.38 overlap score). The overlap scores for the time and population observables do not agree so well, with values zero or near zero. Nevertheless,

this divergence in the population evolution is not critical for the validity of the TD-BA model. As shown in Figure 8, the TD-BA dynamics follow the same qualitative trends in the **h**- and **p**-adjusted datasets. However, the initial TD-BA evolution is slightly slower, 11 ± 1 fs against 9 ± 1 fs of other datasets. Moreover, the TD-BA S₁ population tends to be higher, especially when compared to the **h**-adjusted dataset. When considering all six observables, the mean overlap score Λ⁽¹,²⁾ between TD-BA and the **h**-adjusted reference dataset is 0.40 (Table 3), not much worse than the 0.59 obtained in the comparison between the **p**-adjusted dataset and the **h**-adjusted dataset.

As a final remark, in Ref. 32, Ibele and Curchod reported surface hopping dynamics results for fulvene, also computed at CAS(6,6) level. Their S₁ population at 20 fs computed with **p**-adjusted DC-FSSH was, like ours, near 0.4 (see Figure 4 of that paper). Nevertheless, with **h**-adjusted DC-FSSH, this observable was reduced to about 0.2, which disagrees with our results. Although we cannot provide a definitive explanation for the divergence, it may be that their trajectory ensemble
<figure><img src="image_15.png" /><figcaption>**Figure 8.** Fulvene S<sub>1</sub> population. Top: time-dependent Baeck-An (TD-BA) and **h**-adjusted dataset; bottom: TD-BA and **p**-adjusted dataset. In each case, the shaded area indicates the 95% confidence interval sampled with bootstrap through 10,000 repetitions of the 200 trajectories of each dataset.</figcaption></figure>

composed of 18 initial conditions repeated 10 times each is too biased to deliver a reliable population.

## VIII. Conclusions

In Ref. 27, Baeck and An showed that the first-order nonadiabatic coupling vector could be approximated by a 1D model, depending only on the energy gap and the second derivative of this gap with respect to the coupling coordinate (see Equation (1)). In this paper, the BA approximation is recast on the time domain, giving rise to the TD-BA approximation (Equation (12)), which can be applied to multidimensional systems. The TD-BA couplings are equivalent to a coupling projection on the nuclear velocities, making them particularly suited to be employed in DC-FSSH simulations. TD-BA couplings immediately enable nonadiabatic dynamics with any electronic method that can provide excitation energies and energy gradients.

We have used the methodological protocols proposed in Ref. 39 to evaluate dynamics based on TD-BA couplings. Extensive simulations for ethylene and fulvene at MCSCF level showed that TD-BA DC-FSSH provides a correct qualitative picture of the dynamics and delivered results in semiquantitative agreement with the reference datasets computed with exact nonadiabatic couplings. The results' quality is statistically not much worse than performing DC-FSSH adjusting velocities after hopping in the momentum direction. Nevertheless, the model fails to adequately describe the nonadiabatic dynamics in regions with strong nonadiabatic couplings but small velocities

(which are themselves challenging regions for surface hopping in general<sup>84</sup>).

In our implementation of the DC-FSSH with TD-BA couplings, we tested two ways of computing second time derivatives of the energy gaps and proposed three conditions to be satisfied to reduce numerical instabilities in the calculations. Although this initial modeling has delivered adequate results, there is a margin for further improvements, especially in detecting artifacts caused by small discontinuities in the potential energy surface.

We have also tested the TD-BA DC-FSSH with TDDFT to check whether this method would allow predicting hoppings to the ground state. In this case, however, the answer was negative due to singlet instabilities. This problem is not restricted to TD-BA and should happen in TDDFT dynamics based on exact couplings too. Dynamics based on TD-BA with TDDFT may still help describe internal conversion between excited states.

Given the uncertainties introduced by the TD-BA approximation, DC-FSSH dynamics with these couplings should be restricted to exploratory dynamics, where high accuracy is not a strong requirement. In this work, we have tested the model for the nonadiabatic dynamics of ethylene and fulvene. Although together these two molecules are representative of many typical excited-state cases, dynamics with TD-BA couplings should be thoroughly tested before being applied to other types of topography and densities, like in superexchange<sup>85</sup>, spin-exchange<sup>86</sup>, solvent-solute charge-transfer<sup>87</sup>, three-states crossing<sup>88</sup>, and dissociative<sup>89</sup> internal conversion.

## Data availability

### Underlying data

Figshare: Fulvene DC-FSSH. https://doi.org/10.6084/m9 figshare.14446998.v1<sup>44</sup>.

This project contains the following underlying data:

- All Newton-X input and output files for dynamics of fulvene (text format) with TD-BA, **h**-adjusted velocities, and **p**-adjusted velocities.

Figshare: TD-BA DC-FSSH Dynamics. https://doi.org/10.6084/m9 figshare.14010311.v1<sup>45</sup>.

This project contains the following underlying data:

- All Newton-X input and output files for dynamics of ethylene (text format) with TD-BA.

Figshare: ethylene-dyn-2021.tgz. https://doi.org/10.6084/m9 figshare.13522856.v2<sup>46</sup>.

This project contains the following underlying data:

- All Newton-X input and output files for dynamics of ethylene (text format) with **h**-adjusted velocities and **p**-adjusted velocities.

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).
References

1. Nelson TR, White AJ, Bjorgaard JA, et al.: *Non-adiabatic Excited-State Molecular Dynamics: Theory and Applications for Modeling Photophysics in Extended Molecular Materials.* *Chem Rev.* 2020; **120**(4): 2215-2287.
*PubMed Abstract* | *Publisher Full Text*

2. Crespo-Otero R, Barbatti M: *Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics*. *Chem Rev.* 2018; **118**(15): 7026-7068.
*PubMed Abstract* | *Publisher Full Text*

3. Malone W, Nebgen B, White A, et al.: NEXMD Software Package for Non-adiabatic Excited State Molecular Dynamics Simulations. *J Chem Theory Comput.* 2020; **16**(9): 5771-5783.
*PubMed Abstract* | *Publisher Full Text*

4. Weingart O, Nenov A, Altoè P, et al.: COBRAMM 2.0 - A software interface for tailoring molecular electronic structure calculations and running nanoscale (QM/MM) simulations. *J Mol Model.* 2018; **24**(9): 271.
*PubMed Abstract* | *Publisher Full Text*

5. Humeniuk A, Mitric R: DFTBaby: A Software Package for Non-adiabatic Molecular Dynamics Simulations Based on Long-Range Corrected Tight-binding TD-DFT(B). *Comput Phys Commun.* 2017; **221**: 174-202.
*Publisher Full Text*

6. Du L, Lan Z: *An On-the-Fly Surface-Hopping Program JADE for Nonadiabatic Molecular Dynamics of Polyatomic Systems: Implementation and Applications*. *J Chem Theory Comput.* 2015; **11**(4): 1360-1374.
*PubMed Abstract* | *Publisher Full Text*

7. Barbatti M, Ruckenbauer M, Plasser F, et al.: Newton-X: a surface-hopping program for nonadiabatic molecular dynamics. *WIREs: Comp Mol Sci.* 2014; **4**(1): 26-33.
*Publisher Full Text*

8. Akimov AV, Prezhdo OV: *The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems*. *J Chem Theory Comput.* 2013; **9**(11): 4959-4972.
*PubMed Abstract* | *Publisher Full Text*

9. Richter M, Marquetand P, González-Vázquez J, et al.: SHARC: Ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings. *J Chem Theory Comput.* 2011; **7**(5): 1253-1258.
*PubMed Abstract* | *Publisher Full Text*

10. Send R, Furche F: *First-order nonadiabatic couplings from time-dependent hybrid density functional response theory: Consistent formalism, implementation, and performance*. *J Chem Phys.* 2010; **132**(4): 044107.
*PubMed Abstract* | *Publisher Full Text*

11. Dallos M, Lischka H, Shepard R, et al.: Analytic Evaluation of Nonadiabatic Coupling Terms at the MR-CI Level. II. Minima on the Crossing Seam: Formaldehyde and the Photodimerization of Ethylene. *J Chem Phys.* 2004; **120**(16): 7330-7339.
*PubMed Abstract* | *Publisher Full Text*

12. Park JW, Shiozaki T: *Analytical Derivative Coupling for Multistate CASPT2 Theory*. *J Chem Theory Comput.* 2017; **13**(6): 2561-2570.
*PubMed Abstract* | *Publisher Full Text*

13. Gonon B, Perveaux A, Gatti F, et al.: On the applicability of a wavefunction-free, energy-based procedure for generating first-order non-adiabatic couplings around conical intersections. *J Chem Phys.* 2017; **147**(11): 114114.
*PubMed Abstract* | *Publisher Full Text*

14. Köppel H, Schubert B: *The concept of regularized diabatic states for a general conical intersection*. *Mol Phys.* 2006; **104**: 1069-1079.
*Publisher Full Text*

15. Yu L, Xu C, Lei Y, et al.: *Trajectory-based nonadiabatic molecular dynamics without calculating nonadiabatic coupling in the avoided crossing case: trans--cis photoisomerization in azobenzene.* *Phys Chem Chem Phys.* 2014; **16**(47): 25883-25895.
*PubMed Abstract* | *Publisher Full Text*

16. Zhu C, Nobusada K, Nakamura H: *New implementation of the trajectory surface hopping method with use of the Zhu-Nakamura theory*. *J Chem Phys.* 2001; **115**(7): 3031-3044.
*Publisher Full Text*

17. Suchan J, Janos J, Slavicek P: *Pragmatic Approach to Photodynamics: Mixed Landau-Zener Surface Hopping with Intersystem Crossing*. *J Chem Theory Comput.* 2020; **16**(9): 5809-5820.
*PubMed Abstract* | *Publisher Full Text*

18. Xie W, Sapunar M, Došlic N, et al.: *Assessing the performance of trajectory surface hopping methods: Ultrafast internal conversion in pyrazine*. *J Chem Phys.* 2019; **150**(15): 154119.
*PubMed Abstract* | *Publisher Full Text*

19. Belyaev AK, Lasser C, Trigila G: *Landau-Zener Type Surface Hopping Algorithms*. *J Chem Phys.* 2014; **140**(22): 224108.
*PubMed Abstract* | *Publisher Full Text*

20. Tully JC: *Molecular-Dynamics with Electronic-Transitions*. *J Chem Phys.* 1990; **93**: 1061-1071.
*Publisher Full Text*

21. Hammes-Schiffer S, Tully JC: *Proton Transfer in Solution: Molecular Dynamics with Quantum Transitions*. *J Chem Phys.* 1994; **101**(6): 4657-4667.
*Publisher Full Text*

22. Plasser F, Ruckenbauer M, Mai S, et al.: Efficient and Flexible Computation of Many-Electron Wave Function Overlaps. *J Chem Theory Comput.* 2016; **12**(3): 1207-1219.
*PubMed Abstract* | *Publisher Full Text* | *Free Full Text*

23. Ryabinkin IG, Nagesh J, Izmaylov AF: *Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods*. *J Phys Chem Lett.* 2015; **6**(21): 4200-4203.
*PubMed Abstract* | *Publisher Full Text*

24. Meek GA, Levine BG: *Evaluation of the Time-Derivative Coupling for Accurate Electronic State Transition Probabilities from Numerical Simulations*. *J Phys Chem Lett.* 2014; **5**(13): 2351-2356.
*PubMed Abstract* | *Publisher Full Text*

25. Granucci G, Persico M, Togniolo A: *Direct Semiclassical Simulation of Photochemical Processes with Semiempirical Wave Functions*. *J Chem Phys.* 2001; **114**(24): 10608-10615.
*Publisher Full Text*

26. Westermayr J, Gastegger M, Marquetand P: *Combining SchNet and SHARC: The SchNarc Machine Learning Approach for Excited-State Dynamics*. *J Phys Chem Lett.* 2020; **11**(10): 3828-3834.
*PubMed Abstract* | *Publisher Full Text* | *Free Full Text*

27. Baeck KK, An H: *Practical approximation of the non-adiabatic coupling terms for same-symmetry interstate crossings by using adiabatic potential energies only*. *J Chem Phys.* 2017; **146**(6): 064107.
*PubMed Abstract* | *Publisher Full Text*

28. An H, Baeck KK: *Practical and reliable approximation of nonadiabatic coupling terms between triplet electronic states using only adiabatic potential energies*. *Chem Phys Lett.* 2018; **696**: 100-105.
*Publisher Full Text*

29. Westermayr J, Faber FA, Christensen AS, et al.: *Neural networks and kernel ridge regression for excited states dynamics of CH<sub>3</sub>NH<sub>2</sub>: From single-state to multi-state representations and multi-property machine learning models*. *Machine Learn Sci Technol.* 2020; **1**(1): 025009.
*Publisher Full Text*

30. Hu D, Xie Y, Li X, et al.: *Inclusion of Machine Learning Kernel Ridge Regression Potential Energy Surfaces in On-the-Fly Nonadiabatic Molecular Dynamics Simulation*. *J Phys Chem Lett.* 2018; **9**(11): 2725-2732.
*PubMed Abstract* | *Publisher Full Text*

31. Dral PO, Barbatti M, Thiel W: *Nonadiabatic Excited-State Dynamics with Machine Learning*. *J Phys Chem Lett.* 2018; **9**(19): 5660-5663.
*PubMed Abstract* | *Publisher Full Text* | *Free Full Text*

32. Ibele LM, Curchod BFE: *A molecular perspective on Tully models for nonadiabatic dynamics*. *Phys Chem Chem Phys.* 2020; **22**(27): 15183-15196.
*PubMed Abstract* | *Publisher Full Text*

33. Gozem S, Melaccio F, Valentini A, et al.: *Shape of Multireference, Equation-of-Motion Coupled-Cluster, and Density Functional Theory Potential Energy Surfaces at a Conical Intersection*. *J Chem Theory Comput.* 2014; **10**(8): 3074-3084.
*PubMed Abstract* | *Publisher Full Text*

34. Abiola TT, Rodrigues NDN, Ho C, et al.: *New Generation UV-A Filters: Understanding Their Photodynamics on a Human Skin Mimic*. *J Phys Chem Lett.* 2021; **12**(1): 337-344.
*PubMed Abstract* | *Publisher Full Text*

35. Serrano-Andrés L, Merchán M, Borin AC: *Adenine and 2-aminopurine: Paradigms of modern theoretical photochemistry*. *Proc Natl Acad Sci U S A.* 2006; **103**(23): 8691-8696.
*PubMed Abstract* | *Publisher Full Text* | *Free Full Text*

36. Blancafort L: *Photochemistry and Photophysics at Extended Seams of Conical Intersection*. *Chemphyschem.* 2014; **15**(15): 3166-3181.
*PubMed Abstract* | *Publisher Full Text*

37. Köppel H, Gronki J, Mahapatra S: *Construction scheme for regularized diabatic states*. *J Chem Phys.* 2001; **115**(6): 2377-2388.
*Publisher Full Text*

38. Granucci G, Persico M: *Critical appraisal of the fewest switches algorithm for surface hopping*. *J Chem Phys.* 2007; **126**(13): 134114.
*PubMed Abstract* | *Publisher Full Text*

39. Barbatti M: *Velocity Adjustment in Surface Hopping: Ethylene as a Case Study of the Maximum Error Caused by Direction Choice*. *J Chem Theory Comput.* 2021.
*PubMed Abstract* | *Publisher Full Text*

40. Ditchfield R, Hehre WJ, Pople JA: *Self-Consistent Molecular-Orbital Methods. IX. an Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules*. *J Chem Phys.* 1971; **54**(2): 724-728.
*Publisher Full Text*

41. Lischka H, Shepard R, Müller T, et al.: *The generality of the GUGA MRCI*
approach in COLUMBUS for treating complex quantum chemistry. *J Chem Phys.* 2020; **152**(13): 134110.

PubMed Abstract | Publisher Full Text

42. Lischka H, Shepard R, Shavitt I, et al.: COLUMBUS, an ab initio electronic structure program, release 7.0. 2017. Reference Source

43. Barbatti M, Granucci G, Ruckenbauer M, et al.: NEWTON-X: A package for Newtonian Dynamics Close to the Crossing Seam (v. 2.2). 2018.

44. do Casal MT, Toldo JM, Pinheiro M, et al.: Fulvene DC-FSSH. figshare. Dataset. 2021. http://www.doi.org/10.6084/m9 figshare.14446998.v1

45. Barbatti M: Ethylene DC-FSSH: TD-BA. figshare. Dataset. 2021. http://www.doi.org/10.6084/m9 figshare.14010311.v1

46. Barbatti M: Ethylene DC-FSSH: Momentum Adjustment. figshare. Dataset. 2021. http://www.doi.org/10.6084/m9 figshare.13522856.v2

47. Yanai T, Tew DP, Handy NC: A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP). *Chem Phys Lett.* 2004; **393**(1-3): 51-57.
Publisher Full Text

48. Chai JD, Head-Gordon M: Long-Range Corrected Hybrid Density Functionals With Damped Atom-Atom Dispersion Corrections. *Phys Chem Chem Phys.* 2008; **10**(44): 6615-6620.
PubMed Abstract | Publisher Full Text

49. Balasubramani SG, Chen GP, Coriani S, et al.: TURBOMOLE: Modular program suite for ab initio quantum-chemical and condensed-matter simulations. *J Chem Phys.* 2020; **152**(18): 184107.
PubMed Abstract | Publisher Full Text | Free Full Text

50. TURBOMOLE V7.5 2020, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since. 2007.

51. Barca GMJ, Bertoni C, Carrington L, et al.: Recent developments in the general atomic and molecular electronic structure system. *J Chem Phys.* 2020; **152**(15): 154102.
PubMed Abstract | Publisher Full Text

52. Dawes R, Jasper AW, Tao C, et al.: Theoretical and Experimental Spectroscopy of the S2 State of CHF and CDF: Dynamically Weighted Multireference Configuration Interaction Calculations for High-Lying Electronic States. *J Phys Chem Lett.* 2010; 1: 641-646.
Publisher Full Text

53. Huix-Rotllant M, Ferré N, Barbatti M: Time-Dependent Density Functional Theory: from the Fundamentals to Nonadiabatic Dynamics. In *Quantum Chemistry and Dynamics of Excited States: Methods and Applications*, González, L Lindh, R., Eds. John Wiley & Sons. 2020; 15-46.
Reference Source

54. Casida ME, Huix-Rotllant M: Progress in Time-Dependent Density-Functional Theory. *Annu Rev Phys Chem.* 2012; **63**: 287-323.
PubMed Abstract | Publisher Full Text

55. Levine BG, Ko C, Quenneville J, et al.: Conical intersections and double excitations in time-dependent density functional theory. *Mol Phys.* 2006; **104**(5-7): 1039-1051.
Publisher Full Text

56. Hirata S, Head-Gordon M: Time-Dependent Density Functional Theory within the Tamm-dancoff Approximation. *Chem Phys Lett.* 1999; **314**(3-4): 291-299.
Publisher Full Text

57. Tapavicza E, Tavernelli I, Rothlisberger U, et al.: Mixed Time-Dependent Density-Functional Theory/Classical Trajectory Surface Hopping Study of Oxirane Photochemistry. *J Chem Phys.* 2008; **129**(12): 124108.
PubMed Abstract | Publisher Full Text

58. Fazzi D, Barbatti M, Thiel W: Modeling Ultrafast Exciton Deactivation in Oligothiophenes via Nonadiabatic Dynamics. *Phys Chem Chem Phys.* 2015; **17**(12): 7789-7799.
PubMed Abstract | Publisher Full Text

59. Dreuw A, Wormit M: The Algebraic Diagrammatic Construction Scheme for the Polarization Propagator for the Calculation of Excited States. *WIREs: Comp Mol Sci.* 2015; 5(1): 82-95.
Publisher Full Text

60. Viel A, Krawczyk RP, Manthe U, et al.: Photoinduced dynamics of ethene in the N, V, and Z valence states: a six-dimensional nonadiabatic quantum dynamics investigation. *J Chem Phys.* 2004; **120**(23): 11000-11010.
PubMed Abstract | Publisher Full Text

61. Zhou W, Mandal A, Huo P: Quasi-Diabatic Scheme for Nonadiabatic On-the-Fly Simulations. *J Phys Chem Lett.* 2019; **10**(22): 7062-7070.
PubMed Abstract | Publisher Full Text

62. Ben-Nun M, Quenneville J, Martínez TJ: Ab Initio Multiple Spawning: Photochemistry from First Principles Quantum Molecular Dynamics. *J Phys Chem A.* 2000; **104**(22): 5161-5175.
Publisher Full Text

63. Mori T, Glover WJ, Schuurman MS, et al.: Role of Rydberg States in the Photochemical Dynamics of Ethylene. *J Phys Chem A.* 2012; **116**(11):

2808-2818.

PubMed Abstract | Publisher Full Text

64. Tao H, Allison TK, Wright TW, et al.: Ultrafast internal conversion in ethylene. *J. The excited state lifetime. J Chem Phys.* 2011; **134**(24): 244306.
PubMed Abstract | Publisher Full Text

65. Yu JK, Bannwarth C, Hohenstein EG, et al.: Ab Initio Nonadiabatic Molecular Dynamics with Hole-Hole Tamm-Dancoff Approximated Density Functional Theory. *J Chem Theory Comput.* 2020; **16**(9): 5499-5511.
PubMed Abstract | Publisher Full Text

66. Saita K, Shalashilin DV: On-the-fly ab initio molecular dynamics with multiconfigurational Ehrenfest method. *J Chem Phys.* 2012; **137**(22): 22A506.
PubMed Abstract | Publisher Full Text

67. Barbatti M, Ruckenbauer M, Lischka H: The photodynamics of ethylene: a surface-hopping study on structural aspects. *J Chem Phys.* 2005; **122**(17): 174307.
PubMed Abstract | Publisher Full Text

68. Barbatti M, Granucci G, Persico M, et al.: Semiempirical molecular dynamics investigation of the excited state lifetime of ethylene. *Chem Phys Lett.* 2005; **401**(1-3): 276-281.
Publisher Full Text

69. Mališ M, Luber S: Trajectory Surface Hopping Nonadiabatic Molecular Dynamics with Kohn-Sham ΔSCF for Condensed-Phase Systems. *J Chem Theory Comput.* 2020; **16**(7): 4071-4086.
PubMed Abstract | Publisher Full Text

70. Pradhan E, Sato K, Akimov AV: Non-adiabatic molecular dynamics with ΔSCF excited states. *J Phys Condens Matter.* 2018; **30**(48): 484002.
PubMed Abstract | Publisher Full Text

71. Hollas D, Šišků L, Hohenstein EG, et al.: Nonadiabatic Ab Initio Molecular Dynamics with the Floating Occupation Molecular Orbital-Complete Active Space Configuration Interaction Method. *J Chem Theory Comput.* 2018; **14**(1): 339-350.
PubMed Abstract | Publisher Full Text

72. Sellner B, Barbatti M, Müller T, et al.: Ultrafast Non-Adiabatic Dynamics of Ethylene Including Rydberg States. *Mol Phys.* 2013; **111**(16-17): 2439-2450.
Publisher Full Text

73. Kosma K, Trushin SA, Fuss W, et al.: Ultrafast dynamics and coherent oscillations in ethylene and ethylene-d4 excited at 162 nm. *J Phys Chem A.* 2008; **112**(33): 7514-7529.

PubMed Abstract | Publisher Full Text

74. Kobayashi T, Horio T, Suzuki T: Ultrafast Deactivation of the ππ*(V) State of Ethylene Studied Using Sub-20 fs Time-Resolved Photoelectron Imaging. *J Phys Chem A.* 2015; **119**(36): 9518-9523.

PubMed Abstract | Publisher Full Text

75. Grohmann T, Deeb O, Leibschner M: Quantum separation of para- and ortho-fulvene with coherent light: The influence of the conical intersection. *Chem Phys.* 2007; **338**(2-3): 252-258.
Publisher Full Text

76. Belz S, Grohmann T, Leibschner M: Quantum dynamical simulations for nuclear spin selective laser control of ortho- and para-fulene. *J Chem Phys.* 2009; **131**(3): 034305.
PubMed Abstract | Publisher Full Text

77. Alfalah S, Belz S, Deeb O, et al.: Photoinduced quantum dynamics of ortho- and para-fulvene: Hindered photoisomerization due to mode selective fast radiationless decay via a conical intersection. *J Chem Phys.* 2009; **130**(12): 124318.
PubMed Abstract | Publisher Full Text

78. Blancafort L, Gatti F, Meyer HD: Quantum dynamics study of fulene double bond photoisomerization: The role of intramolecular vibrational energy redistribution and excitation energy. *J Chem Phys.* 2011; **135**(13): 134303.
PubMed Abstract | Publisher Full Text

79. Ruiz-Barragan S, Blancafort L: Photophysics of fulene under the non-resonant stark effect. Shaping the conical intersection seam. *Faraday Discuss.* 2013; **163**: 497-512.
PubMed Abstract | Publisher Full Text

80. Mendive-Tapia D, Lasorne B, Worth GA, et al.: Controlling the mechanism of fulene S(1)/S(0) decay: switching off the stepwise population transfer. *Phys Chem Chem Phys.* 2010; **12**(48): 15725-15733.
PubMed Abstract | Publisher Full Text

81. Mendive-Tapia D, Lasorne B, Worth GA, et al.: Towards converging non-adiabatic direct dynamics calculations using frozen-width variational Gaussian product basis functions. *J Chem Phys.* 2012; **137**(22): 22A548.
PubMed Abstract | Publisher Full Text

82. Bearpark MJ, Bernardi F, Olivucci M, et al.: Can Fulvene S, Decay Be Controlled? A CASSIM Study with MMVB Dynamics. *J Am Chem Soc.* 1996; **118**(22): 5254-5260.
Publisher Full Text

83. Bearpark MJ, Blancafort L, Paterson MJ: Mapping the intersection space of the ground and first excited states of fulvene. *Mol Phys.* 2006; **104**(5-7): 1033-1038.
Publisher Full Text

84. Subotnik JE, Ouyang W, Landry BR: Can We Derive Tully's Surface-Hopping
**Algorithm from the Semiclassical Quantum Liouville Equation? Almost, but Only with Decoherence.** *J Chem Phys.* 2013; **139**(21): 214107.
**PubMed Abstract** | **Publisher Full Text**

85. Wang L, Trivedi D, Prezhdo OV: **Global Flux Surface Hopping Approach for Mixed Quantum-Classical Dynamics.** *J Chem Theory Comput.* 2014; **10**(9): 3598–3605.
**PubMed Abstract** | **Publisher Full Text**

86. Bai S, Barbatti M: **Mechanism of Spin-Exchange Internal Conversion: Practical Proxies for Diabatic and Nonadiabatic Couplings.** *J Chem Theory Comput.* 2019; **15**(3): 1503–1513.
**PubMed Abstract** | **Publisher Full Text**

87. Barbatti M: Photorelaxation Induced by Water-Chromophore Electron Transfer. *J Am Chem Soc.* 2014; **136**(29): 10246–10249.
**PubMed Abstract** | **Publisher Full Text**

88. Kistler KA, Matsika S: **Three-state conical intersections in cytosine and pyrimidinone bases.** *J Chem Phys.* 2008; **128**(21): 215102.
**PubMed Abstract** | **Publisher Full Text**

89. Sobolewski AL, Domcke W, Dedonder-Lardeux C, et al.: **Excited-State Hydrogen Detachment and Hydrogen Transfer Driven by Repulsive $^1\pi\sigma^*$ States: A New Paradigm for Nonradiative Decay in Aromatic Biomolecules.** *Phys Chem Chem Phys.* 2002; **4**(7): 1093–1100.
**Publisher Full Text**
# Open Peer Review

Current Peer Review Status: ✓ ✓ ✓

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

## Version 1

### Reviewer Report 02 June 2021

https://doi.org/10.21956/openreseurope.14695.r26914

© 2021 **Blancafort L.** This is an open access peer review report distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_19.png" /></figure>

### Lluís Blancafort

<figure><img src="image_20.png" /></figure>

Department of Chemistry and Institute of Computational Chemistry and Catalysis. Faculty of Sciences, University of Girona, Girona, Spain

This manuscript introduces a new approximation to calculate non-adiabatic couplings for trajectory surface hopping, based on an expression introduced previously by Baeck and An. The development of new schemes for coupling calculation is important to make excited-state simulations more efficient. The authors provide a detailed description of their new approximation that ensures that it can be fully reproduced. The same can be said about the applications (fulvene and ethylene excited state trajectories). Therefore, the work is fully suitable for indexing in Open Research Europe.

As a minor revision, I would ask the authors to comment on the validity of the assumption that the only non-null energy Hessian component is the one in the Q direction, where Q is the coupling coordinate (see Section III). This assumption seems key to extend the BA's approach to multidimensional systems. However in polyatomic molecules the opposite often occurs, i.e. the S0 and S1 surface have different curvature along several vibrational coordinates. This has been studied in detail e.g. for fulvene¹, and large differences in the S0 and S1 Hessian are found along symmetry-preserving coordinates where the coupling is zero. Will the new approach not overestimate the coupling in such cases?

### References

1. Paterson MJ, Bearpark MJ, Robb MA, Blancafort L: The curvature of the conical intersection seam: an approximate second-order analysis. *J Chem Phys*. 2004; **121** (23): 11562-71 [PubMed Abstract](https://pubmed.com/11562-71) | [Publisher Full Text](https://www pubmed com/11562-71)

**Is the rationale for developing the new method (or application) clearly explained?**

Yes

**Is the description of the method technically sound?**

Yes
**Are sufficient details provided to allow replication of the method development and its use by others?**

Yes

**If any results are presented, are all the source data underlying the results available to ensure full reproducibility?**

Yes

**Are the conclusions about the method and its performance adequately supported by the findings presented in the article?**

Yes

**Competing Interests:** No competing interests were disclosed.

**Reviewer Expertise:** Computational chemistry.

**I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.**

Reviewer Report 01 June 2021

https://doi.org/10.21956/openreseurope.14695.r26913

© 2021 Westermayr J. This is an open access peer review report distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_21.png" /></figure>

Julia Westermayr

<figure><img src="image_22.png" /></figure>

Computational Surface Chemistry Group, Department of Chemistry, University of Warwick, Coventry, UK

The article "Fewest switches surface hopping with Baeck-An couplings" by Mariana T. do Casal et al. represents a detailed analysis of the Baeck-An couplings within a time-dependent approximation (TD-BA) that can be used in fewest-switches surface hopping molecular dynamics simulations. This approach is especially powerful as it can enable surface hopping dynamics using only predicted energies and forces. Their analysis is based on the ultrafast dynamics of ethylene and fulvene at MCSCF level, two molecules recently considered by Ibele et al. as molecular Tully models for benchmarking different nonadiabatic dynamics methods. The authors compare their results with nonadiabatic dynamics with "exact" nonadiabatic coupling vectors, i.e., computed with the reference method, in combination with velocities adjusted in momentum and nonadiabatic coupling direction. Their results show that the TD-BA approximation leads to qualitatively correct trends in the dynamics and a semiquantitative agreement with respect to the reference dynamics. The authors further tested the method for fulvene and thiophene at the TDDFT level of theory. However, the use of the TD-BA model is restricted to excited state hoppings. Hoppings from an excited state to the ground state remain problematic, but this problem is also experienced in surface hopping dynamics with "exact" nonadiabatic coupling vectors. Their work clearly
demonstrates that the TD-BA approximation can become a comparative method for studying the photodynamics of molecular systems when an accurate reference method is available and clearly highlights the discrepancies and uncertainties of the different flavours of surface hopping dynamics.

Overall, the authors have presented their manuscript in a clearly structured manner and I enjoyed reading the article. The manuscript is well written, it is easy to follow how conclusions are drawn and explanations are provided when needed. I especially think that the authors have done a fantastic job of laying out the technical details in Section VI A. This section explains how to obtain the second time derivative of the energy gap and what approaches are tested to address the discontinuities in the potential energy landscapes introduced by the reference method.

I have a few minor suggestions that can further improve the manuscript:

1. To analyse the differences between the methods, the authors use an (average) overlap score that was also used in Ref. 39, in which the corresponding author compared the surface hopping dynamics of ethylene with different types of velocity adjustment. Twice in the text, the reader is referred to Ref. 39. I believe that a brief description of the method (in somewhat more detail than in the second paragraph on page 6) would benefit readers.

2. The authors describe their efforts to simulate the dynamics of fulvene using TDDFT and TDA-DFT. Due to problems in the reference simulations, the authors tested their method for the surface hopping dynamics of thiophene. Besides this reason for studying thiophene, it would be interesting to know what other reasons led to the choice of thiophene. A brief description of its dynamics and why it is particularly suitable for the purpose of this study would be beneficial.

3. Ibele et al. (ref. 32) proposed three Tully models, i.e., ethylene, fulvene, and DMABN, to benchmark surface hopping methods. Is there any reason why the authors have chosen ethylene and fulvene and did not consider DMABN? I do not expect the authors to perform their analysis for the third molecule, but I believe that a sentence or two in which the authors explain their choice might be helpful. I would also appreciate it if the authors could comment on the expected performance of the TD-BA approximation in studying the dynamics of DMABN.

4. Page 6, first paragraph: I believe there is a typo and "εω" should read "δω".

5. Page 8, Fig. 4 (d): The authors state that "the TD-BA model shows spurious spikes that even the cleaning conditions discussed above could not get rid of." However, looking at Fig. 4(a) suggests that the spikes between 60 and 80 fs could be eliminated with a small enough value of δε as the energy gap between S1 and S2 seems to be in the range of 5 eV.

6. If the authors wish, they could slightly improve Fig. 7, as I think this graph could benefit from molecular structures of fulvene.

Besides these minor comments, I do not have any further suggestions and believe that the article meets the criteria for indexing in Open Research Europe.

**Is the rationale for developing the new method (or application) clearly explained?**
Yes

**Is the description of the method technically sound?**

Yes

**Are sufficient details provided to allow replication of the method development and its use by others?**

Yes

**If any results are presented, are all the source data underlying the results available to ensure full reproducibility?**

Yes

**Are the conclusions about the method and its performance adequately supported by the findings presented in the article?**

Yes

**Competing Interests:** No competing interests were disclosed.

**Reviewer Expertise:** Nonadiabatic molecular dynamics, surface hopping, nonadiabatic coupling approximation, machine learning.

**I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.**

Reviewer Report 25 May 2021

https://doi.org/10.21956/openreseurope.14695.r26916

© 2021 Fernadez-Alberti S. This is an open access peer review report distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_23.png" /></figure>

Sebastian Fernadez-Alberti

Universidad Nacional de Quilmes, Bernal, Argentina

The authors statistically analyzed the different approximations by which the Baeck-An couplings can be implemented is FSSH algorithms. They have tested the method using ethylene and fulvene molecules. They provide approximations given in eq. 12 and also the calculation of FSSH probability(eq. 13). Additional TD-BA DC-FSSH calculations were also done with linear-response time-dependent density functional (TDDFT). Besides, the authors explain how the second time derivative in Equation (12) can be computed numerically or analytically using eq. 15-17. The article property explains the different approximation, achieved accuracies and benchmark tests. They proved that TD-BA couplings are equivalent to a coupling projection on the nuclear velocities, making them particularly suited to be employed in DC-FSSH simulations. Extensive simulations for ethylene and fulvene at MCSCF level showed that TD-BA combined with DC-FSSH provide a correct
qualitative picture of the dynamics and delivered results in semiquantitative agreement with the reference datasets computed with exact nonadiabatic couplings. The results' quality is statistically not much worse than performing DC-FSSH adjusting velocities after hopping in the momentum direction. Nevertheless, the model fails to adequately describe the nonadiabatic dynamics in regions with strong nonadiabatic couplings but small velocities. Finally, the authors test two different ways of computing second time derivatives of the energy gaps and proposed three conditions to be satisfied to reduce numerical instabilities in the calculations.

Although this initial modeling has delivered adequate results, the authors claim that there is a margin for further improvements, especially in detecting artifacts caused by small discontinuities in the potential energy surface. The article is well written. I encourage the authors to describe the further improvements they have made and, in case they have performed preliminary tests, perhaps they can describe their findings and, in this way, help other authors future attempts of implementations. Despite that I think the article merits indexing at Open Research Europe.

**Is the rationale for developing the new method (or application) clearly explained?**

Yes

**Is the description of the method technically sound?**

Yes

**Are sufficient details provided to allow replication of the method development and its use by others?**

Yes

**If any results are presented, are all the source data underlying the results available to ensure full reproducibility?**

Yes

**Are the conclusions about the method and its performance adequately supported by the findings presented in the article?**

Yes

**Competing Interests:** No competing interests were disclosed.

**I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.**

Author Response 01 Jun 2021

**Mario Barbatti**, Aix Marseille University, CNRS, ICR, Marseille, France

We thank the reviewer for the report on our paper. Concerning the comment: "*I encourage the authors to describe the further improvements they have made and, in case they have performed preliminary tests, perhaps they can describe their findings and, in this way, help other authors future attempts of implementations.*" At this moment, we still do not have any advances to report. Nevertheless, an advantage of publishing in a platform like Open
Research Europe is that it is possible to update the paper to reflect new developments. We intend to report here any further methodological advances, as suggested by the reviewer.

**Competing Interests:** No competing interests were disclosed.