# Investigating the photodynamics of trans-azobenzene with coupled trajectories

Carlotta Pieroni, Eduarda Sangiogo Gil, Lea Ibele, Maurizio Persico, Giovanni Granucci, Federica Agostini

► To cite this version:

Carlotta Pieroni, Eduarda Sangiogo Gil, Lea Ibele, Maurizio Persico, Giovanni Granucci, et al. Investigating the photodynamics of trans-azobenzene with coupled trajectories. Journal of Chemical Theory and Computation, 2024, 20 (2), pp.580-596. 10.1021/acs.jctc.3c00978. hal-04443966

HAL Id: hal-04443966

https://universite-paris-saclay.hal.science/hal-04443966v1

Submitted on 7 Feb 2024

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
Investigating the photodynamics of *trans*-azobenzene with coupled
trajectories

Carlotta Pieroni¹,²*, Eduarda Sangiogo Gil²*,†, Lea M. Ibele², Maurizio Persico¹, Giovanni Granucci¹ and Federica Agostini²‡

¹University of Pisa, Dipartimento di Chimica e Chimica Industriale, via G. Moruzzi 13,
56124, Pisa, Italy

²Université Paris-Saclay, CNRS, Institut de Chimie Physique UMR8000, 91405, Orsay,
France

September 6, 2023

Abstract

In this work, we present the first implementation of the coupled-trajectory Tully surface hopping (CT-TSH) suitable for applications to molecular systems. We combine CT-TSH with the semiempirical Floating Occupation Molecular Orbitals-Configuration Interaction (FOMO-CI) electronic structure method to investigate the photoisomerization dynamics of *trans*-azobenzene. Our study shows that CT-TSH can capture correctly decoherence effects in this system, yielding consistent electronic and nuclear dynamics in agreement with (standard) decoherence-corrected TSH. Specifically, CT-TSH is derived from the exact factorization and the electronic coefficients' evolution is directly influenced by the coupling of trajectories, resulting in the improvement of internal consistency if compared to standard TSH.

# 1 Introduction

The time evolution of a (non-relativistic) molecular system is determined by the time-dependent Schrödinger equation (TDSE). Unfortunately, solving the TDSE and simulating quantum dynamics in molecules, especially when electronically excited states are involved, remains a fundamentally difficult problem, due to the non-local nature of quantum mechanics and the exponential scaling of the computational effort with the number of degrees of freedom. One popular and effective approximation to circumvent these issues is to

<footnote>*These authors contributed equally.</footnote>

<footnote>†eduarda.sangiogo-gil@universite-paris-saclay.fr</footnote>

<footnote>‡federica.agostini@universite-paris-saclay.fr</footnote>
describe the nuclear motion with classical mechanics instead of quantum mechanics [1, 2, 3, 4, 5, 6, 7, 8]. However, as a consequence of such a classical approximation, some purely quantum effects are lost or improperly described. Quantum decoherence is the main example, as well as quantum interference [9, 10, 11, 12, 7]. In this respect, the quantum-classical numerical approaches derived from the exact factorization seem to show an improvement in the description of quantum decoherence, compared to, for example, the Tully surface hopping (TSH) method (one of most widely used quantum-classical methods in the community) [13]. Exact-factorization-based approaches allow to recover some quantum effects thanks to the coupling among the trajectories [14, 15, 16]. In particular, the recently developed coupled-trajectory Tully surface hopping (CT-TSH) method combines the basic idea of a surface-hopping procedure, i.e., classical adiabatic nuclear evolution using stochastic hops, with the “decoherence-corrected” electronic quantum dynamics, where the additional term accounting for decoherence effects is rigorously derived from the exact factorization [17].

TSH has a systematic problem due to the disconnect between how electrons and nuclei evolve, an issue commonly referred to as “overcoherence”: throughout time, the nuclei evolve on a single Born-Oppenheimer (BO) potential energy surface (PES) but can, at any time, hop to another PES according to a stochastic algorithm, while the electronic evolution remains in a coherent superposition of BO states. To demonstrate the problem of overcoherence, let’s examine the scheme shown in Fig. 1. In this scheme, a trajectory moves from left to right in one dimension and passes through two consecutive avoided crossings of BO PESs. The trajectory starts on the excited state $S_1$, with 100% of the electronic populations corresponding to that state ($t = 0$). As the trajectory passes through the first avoided crossing ($t = t'$), 70% of its population is transferred to the lower state $S_0$ and a surface hop occurs. From that point on, the trajectory follows the $S_0$ PES, and hence the trajectory is able to move further to the right. Without the hop, the trajectory would have been reflected by a classically forbidden barrier in the upper state. The central problem is that the 30% of the population remaining in $S_1$ also follows the gradient of the lower state and moves to the right, although this should be classically forbidden. When the trajectory reaches the second crossing point ($t = t''$), unphysical interference occurs between the 30% of population on the upper state and the 70% on the lower state. Such interference might strongly affect the transferred population depending on the relative complex phase between the two states. As the upper amplitudes should not have reached this crossing at all, no interference should occur. To overcome this inconsistency, several decoherence corrections have been proposed, which, like the TSH procedure itself, are somewhat *ad hoc*, even if physically motivated [18, 19, 20, 21, 22, 23, 24]. Usually, those corrections have the effect of reducing the difference between the quantum electronic populations, the $|C|^2$ of Fig. 1, and the classical populations, obtained by running many trajectories and determining the fraction associated to each state. In the following, these two populations will be indicated as $P(t)$ and $F(t)$, respectively. Instead, in the CT-TSH algorithm, the electronic evolution is directly affected by the coupling of trajectories, which naturally arises in the classical limit of the nuclear dynamics in the exact factorization. As a result, decoherence naturally is included, thereby improving internal consistency without requiring additional *ad hoc* corrections. The coupling among the trajectories
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of overcoherence in a system with a double crossing between the $S_0$ (red) and $S_1$ (blue) BO PESs. The trajectory starts on the upper state, undergoes a surface hop to the lower state, but the remaining amplitude on the upper state also moves to the right, violating classical rules. At the second crossing, interference occurs between the upper and lower populations.</figcaption></figure>

requires, however, a larger computational effort if compared to independent-trajectory TSH, which is a trivially parallelizable algorithm.

Until now, applications of CT-TSH have been limited to model systems of one or two dimensions including only two electronic states [17]. In this study, we aim to demonstrate its efficiency and reliability for the first time on a molecular system. We investigate, in full dimensionality, the photoisomerization of the *trans*-azobenzene (TAB) molecule combining the semiempirical FOMO-CI electronic structure [25] method with the CT-TSH algorithm. We choose to use the semiempirical Floating Occupation Molecular Orbitals-Configuration Interaction (FOMO-CI) method because it offers high computational efficiency and reasonable accuracy (when properly parametrized), making it ideal for easily computing multiple trajectories for medium to large molecular systems over time scales of up to few picoseconds [26, 27, 28, 29].

This paper is organized as follows. First, in Sec. 2.1, we provide a brief overview of the exact factorization method. Next, in Sec. 2.2, we describe the CT-TSH algorithm. We then present the implementation details of CT-TSH for multidimensional systems in Sec. 3, followed by the computational details in Sec. 4. In Sec. 5, we present the results of our study of the photoisomerization dynamics of TAB upon $n\pi^*$ and $\pi\pi^*$ excitation. We conclude our paper with a summary of our findings in Sec. 6.
# 2 Theory and methods

## 2.1 Exact factorization

The interaction between electrons and nuclei in a molecule can be described using the non-relativistic molecular Hamiltonian

$$ \hat{\mathcal{H}} = \sum_{\nu=1}^{N_n} \frac{-\hbar^2}{2M_{\nu}} \nabla_{\nu}^2 + \hat{\mathcal{H}}_{\text{BO}}(\mathbf{r}, \mathbf{R}) \quad (1) $$

where the first term on the right-hand side (RHS) is the nuclear kinetic energy operator expressed in Cartesian coordinates and $\hat{\mathcal{H}}_{\text{BO}}(\mathbf{r}, \mathbf{R})$ is the electronic (or BO) Hamiltonian which contains the electronic kinetic energy and all the interactions. Here, we assigned an index $\nu$ to label the $N_n$ nuclei, while indicating their individual masses as $M_{\nu}$. The symbols $\mathbf{R}$ and $\mathbf{r}$ indicate the sets of $3N_n$ and $3N_e$ nuclear and electronic coordinates, respectively. The time evolution of the electron-nuclear system is dictated by the TDSE

$$ i\hbar \frac{\partial}{\partial t} \Psi(\mathbf{r}, \mathbf{R}, t) = \hat{\mathcal{H}} \Psi(\mathbf{r}, \mathbf{R}, t) \quad (2) $$

whose solution yields the time-dependent molecular wavefunction $\Psi(\mathbf{r}, \mathbf{R}, t)$.

In the exact factorization framework [14, 30, 31], the molecular wavefunction is factored as a single product of a nuclear wavefunction and an electronic conditional factor,

$$ \Psi(\mathbf{r}, \mathbf{R}, t) = \chi(\mathbf{R}, t)\phi(\mathbf{r}, t; \mathbf{R}), \quad (3) $$

where the nuclear wavefunction $\chi(\mathbf{R}, t)$ evolves according to the nuclear TDSE

$$ i\hbar \frac{\partial \chi(\mathbf{R}, t)}{\partial t} = \left[ \sum_{\nu=1}^{N_n} \frac{[-i\hbar\nabla_{\nu} + \mathbf{A}_{\nu}(\mathbf{R}, t)]^2}{2M_{\nu}} + \epsilon(\mathbf{R}, t) \right] \chi(\mathbf{R}, t), \quad (4) $$

and the electronic conditional factor $\phi(\mathbf{r}, t; \mathbf{R})$, with a parametric dependence on $\mathbf{R}$, evolves according to the electronic equation

$$ i\hbar \frac{\partial \phi(\mathbf{r}, t; \mathbf{R})}{\partial t} = [\hat{\mathcal{H}}_{\text{BO}}(\mathbf{r}; \mathbf{R}) + \hat{U}_{\text{en}}[\phi, \chi] - \epsilon(\mathbf{R}, t)] \phi(\mathbf{r}, t; \mathbf{R}). \quad (5) $$

The coupled nuclear and electronic equations of motion (eqs. 4 and 5) describe the dynamical interaction between electrons and nuclei beyond the adiabatic regime. This interaction is mediated by the time-dependent vector potential (TDVP) $\mathbf{A}_{\nu}(\mathbf{R}, t)$, the time-dependent potential energy surface (TDPES) $\epsilon(\mathbf{R}, t)$, and the electron-nuclear coupling operator $\hat{U}_{\text{en}}[\phi, \chi]$. These quantities are defined as:

$$ \mathbf{A}_{\nu}(\mathbf{R}, t) = \langle \phi(t; \mathbf{R}) | -i\hbar\nabla_{\nu} | \phi(t; \mathbf{R}) \rangle_{\mathbf{r}} \quad (6) $$

$$ \epsilon(\mathbf{R}, t) = \left\langle \phi(t; \mathbf{R}) \left| \hat{\mathcal{H}}_{\text{BO}} + \hat{U}_{\text{en}}[\phi, \chi] - i\hbar \frac{\partial}{\partial t} \right| \phi(t; \mathbf{R}) \right\rangle_{\mathbf{r}} \quad (7) $$

and

$$ \hat{U}_{\text{en}}[\phi, \chi] = \sum_{\nu=1}^{N_n} \frac{1}{M_{\nu}} \left[ \frac{[-i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R}, t)]^2}{2} + \left( \frac{-i\hbar\nabla_{\nu}\chi(\mathbf{R}, t)}{\chi(\mathbf{R}, t)} + \mathbf{A}_{\nu}(\mathbf{R}, t) \right) (-i\hbar\nabla_{\nu} - \mathbf{A}_{\nu}(\mathbf{R}, t)) \right] \quad (8) $$
where the symbol $\langle \cdot \rangle_r$ stands for integration over electronic positions.

The ambiguity of the product form of the molecular wavefunction in Eq. 3 is partially eliminated by imposing a partial normalization condition $1 = \langle \phi(t; \mathbf{R}) | \phi(t; \mathbf{R}) \rangle_r \forall \mathbf{R}, t$, which allow us to identify the nuclear marginal probability density as $|\chi(\mathbf{R}, t)|^2 = \langle \Psi(\mathbf{R}, t) | \Psi(\mathbf{R}, t) \rangle_r$ (note that $\Psi$ depends on $\mathbf{r}$, but we use the notation such that the variable that is integrated out does not appear in the bra-ket). As a result, the nuclear and electronic wavefunctions are uniquely defined up to a phase factor $e^{(i/\hbar)\theta(\mathbf{R},t)}$, with $\theta(\mathbf{R}, t)$ being a real function.

## 2.2 Coupled-trajectory approaches

Some of the trajectory-based approaches derived from the exact factorization framework allow propagating coupled electron-nuclear dynamics in the presence of nonadiabatic effects [15, 32, 33, 34, 35, 36] by approximating the nuclear dynamics described by Eq. 4 using an ensemble of coupled trajectories. When the positions and momenta of these trajectories evolve according to Hamilton equations of motion driven by the TDVP and TDPES, the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm is obtained [15, 32, 33]. CT-MQC requires the calculation of the non-adiabatic coupling (NAC) vectors at all times to evaluate the nuclear forces. This is clearly a computationally expensive procedure for medium to large-size molecules. Furthermore, CT-MQC does not guarantee energy conservation either along a single trajectory or over the ensemble [37]. To circumvent these issues, CT-TSH has been introduced: nuclear forces are purely adiabatic, i.e., are obtained as the gradient of a BO PES, and energy conservation is imposed after each hop as in standard TSH along each trajectory. It should be noted that Min and coworkers have proposed alternative trajectory-based procedures derived from the exact factorization but based on independent or auxiliary trajectories [38, 24, 39, 40].

In both CT-MQC and CT-TSH, the electronic wavefunction is expanded in the adiabatic basis, which consists of eigenstates of the electronic Hamiltonian $\hat{H}_{EC}(\mathbf{r}, \mathbf{R})$. The time evolution equations for the expansion coefficients are derived from Eq. 5 and are solved in conjunction with the nuclear Hamilton equations.

From now on, the term “trajectory” will denote a set of $3N_n$ nuclear coordinates that evolve in time, and all quantities in the quantum formulation that depend on $\mathbf{R}$ will depend on the trajectory $\mathbf{R}^{\alpha}(t)$, with $\alpha$ indexing the trajectories. When the electronic wavefunction $\phi(\mathbf{r}, t; \mathbf{R}^{\alpha}(t))$ is expressed as a linear combination of adiabatic states, $\phi^{(m)}(\mathbf{r}, t; \mathbf{R}^{\alpha}(t))$, their parametric dependence on the nuclear position/trajectory induces an implicit time dependence. The expansion coefficients in

$$ \phi(\mathbf{r}, t; \mathbf{R}^{\alpha}(t)) = \sum_{m} C_m(\mathbf{R}^{\alpha}(t), t) \phi^{(m)}(\mathbf{r}, t; \mathbf{R}^{\alpha}(t)) \quad (9) $$

depend on the nuclear trajectory as well. In the following, the dependence on $\mathbf{R}^{\alpha}(t)$ will be indicated only via the index $\alpha$ for simplicity.
In CT-TSH, the trajectories are propagated adiabatically with the force

$$F_{\nu}^{\alpha}(t) = -\nabla_{\nu} E_{*}^{\alpha} \qquad (10)$$

where $E_{*}^{\alpha}$ is the value of the adiabatic energy at the position $R^{\alpha}(t)$ on the active electronic state *, which is selected stochastically at each time step according to the fewest-switches probability [13].

The electronic equation 5 yields the evolution equations for the expansion coefficients $C_{m}^{\alpha}(t)$,

$$\dot{C}_{m}^{\alpha}(t) = \dot{C}_{m,tsh}^{\alpha}(t) + \dot{C}_{m,ct}^{\alpha}(t) \qquad (11)$$

where the first term is the same as in the standard TSH method, i.e.,

$$\dot{C}_{m,tsh}^{\alpha}(t) = -\frac{i}{\hbar} E_{m}^{\alpha} C_{m}^{\alpha}(t) - \sum_{l} \sum_{\nu=1}^{Nn} \mathbf{v}_{\nu}^{\alpha}(t) \cdot \mathbf{d}_{\nu,ml}^{\alpha} C_{l}^{\alpha}(t) \qquad (12)$$

with $\mathbf{v}_{\nu}^{\alpha}(t)$ being the nuclear velocity and $\mathbf{d}_{\nu,ml}^{\alpha} = \langle \varphi^{(m)}(R^{\alpha}(t)) | \nabla_{\nu} \varphi^{(l)}(R^{\alpha}(t)) \rangle$ the NAC vectors. The second term in Eq. 11

$$\dot{C}_{m,ct}^{\alpha}(t) = \sum_{\nu=1}^{Nn} \frac{\mathcal{P}_{\nu}^{\alpha}(t)}{\hbar M_{\nu}} \cdot (\mathbf{f}_{\nu,m}^{\alpha} - \mathbf{A}_{\nu}^{\alpha}(t)) C_{m}(t) \qquad (13)$$

depends on the quantum momentum $\mathcal{P}_{\nu}^{\alpha}(t) = \frac{-\hbar \nabla_{\nu} |\chi(R^{\alpha}(t),t)|^2}{2|\chi(R^{\alpha}(t),t)|^2}$, on the accumulated adiabatic force $\mathbf{f}_{\nu,m}^{\alpha} = \int_{0}^{t} (-\nabla_{\nu} E_{m}^{\alpha}) d\tau$ and on the TDVP. In CT-TSH, the TDVP is approximated as

$$\mathbf{A}_{\nu}^{\alpha}(t) = \hbar \sum_{m,l} \Im[C_{m}^{*\alpha}(t)C_{l}^{\alpha}(t)] \mathbf{d}_{\nu,ml}^{\alpha} + \sum_{m} |C_{m}^{\alpha}(t)|^2 \mathbf{f}_{\nu,m}^{\alpha} \approx \sum_{m} |C_{m}^{\alpha}(t)|^2 \mathbf{f}_{\nu,m}^{\alpha}. \qquad (14)$$

The contribution containing the NAC vectors can be neglected since $\mathbf{d}_{\nu,ml}^{\alpha}$ is usually localized in space, and thus, it is expected to be negligible in comparison to $\mathbf{f}_{\nu,m}^{\alpha}$, which is accumulated in time. This approximation might not be accurate if the trajectories spend a long time in a strong coupling region as the first term on the RHS of Eq. 14 becomes as important, or even larger, than the second term. Relaxing such an approximation in Eq. 14 requires to derive a new implementation for the calculation of the ct term in Eq. 13, which is beyond the scope of this work (but currently under investigation).

Note that the expression of the quantum momentum depends on the changes in the nuclear density across space. In CT-TSH, a semiclassical representation of the nuclear wavefunction is utilized, whereby the nuclear density is reconstructed as a sum of Gaussians centered at the positions of the $N_T$ trajectories. At the position $R^{\alpha}(t)$ of the trajectory $\alpha$, the value of the nuclear density is approximated as

$$|\chi(R^{\alpha}(t), t)|^2 = \frac{1}{N_T} \sum_{\beta=1}^{N_T} G^{\alpha\beta} \qquad (15)$$

where $G^{\alpha\beta}$ is a shorthand for

$$G(R^{\alpha}(t) - R^{\beta}(t)) = \prod_{i=1}^{N_C} \sqrt{\frac{1}{\pi \sigma_i^2}} \exp \left[ \frac{-(R_i^{\alpha}(t) - R_i^{\beta}(t))^2}{\sigma_i^2} \right]. \qquad (16)$$

The symbol $N_C$ stands for the number of degrees of freedom. Further details regarding the choice of the Gaussian widths $\sigma_i$ can be found in Sec. 3.
Eq. 15 demonstrates that in order to calculate the nuclear density, it is necessary to have information about the positions of all trajectories at a given time *t*. Consequently, the trajectories cannot be evolved independently: from the practical point of view, this is the main drawback of CT-TSH if compared to standard TSH.

The CT-TSH method can be an effective alternative to CT-MQC as it employs the TSH scheme to govern nuclear dynamics while still retaining complete coupling among trajectories in the electronic dynamics. In this way, each trajectory is propagated using a purely adiabatic force, while the electronic coefficients are evolved according to CT-MQC. Therefore, CT-TSH trajectories remain coupled through the quantum momentum. The aim of this work is to show that the introduction of the coupled trajectory term of Eq. 13 (i.e., the quantum momentum) allows to correct the overcoherence error of standard TSH [41].

### 3 Implementation of CT-TSH for molecular systems

In this study, we employ the CT-TSH approach for the first time in a molecular application in full dimensionality, to investigate its performance in addressing the challenging issue of overcoherence. Electronic energies, couplings and forces required to propagate the electronic wavefunction and the nuclear trajectories are calculated using a semiempirical configuration interaction (CI) approach based on molecular orbitals (MO) obtained by a self-consistent field (SCF) with floating occupation numbers (FOMO-CI) electronic structure method [25]. Semiempirical methods are known to be significantly faster than *ab initio* ones, which makes them suitable for direct dynamics. With a proper parameterization, these models can reproduce essential features of the excited PESs. Additionally, these models have been validated for several classes of compounds and environments [42, 43, 29].

We adopt the local diabatization (LD) algorithm for the electronic wavefunction propagation, which is numerically stable in strong coupling regions and exempt from the so-called “trivial crossing” problem [25, 2, 44, 45, 46, 47]. Therefore, Eq. 12 is reformulated in the LD formalism as shown in Appendix A. The fourth order Runge-Kutta (RK4) integrator was employed for the propagation of the electronic coefficients, since the presence of the term $\dot{C}_{m,ct}^{\alpha}$ in Eq. 11 prevents us to use the propagator employed in the LD scheme (see Appendix A).

In a nutshell, the CT-TSH algorithm is a variation of the TSH method that includes an extra term, $\dot{C}_{m,ct}^{\alpha}$, in the electronic coefficients' propagations. This term depends on the accumulated adiabatic forces and quantum momentum. The latter is computed from the nuclear density, which is a fundamental quantum mechanical quantity that cannot be exactly calculated in CT-TSH. However, it can be reconstructed using a sum of Gaussians centered at the positions of the nuclear trajectories. The accuracy of the approximation of the nuclear density is crucial for a proper description of the quantum momentum and, consequently, for the propagation of electronic coefficients. It is worth noting that reconstructing the nuclear density for a polyatomic molecule is much more challenging than for a few-dimensional system (as the ones previously
studied with CT-TSH). Hence, different flavours of CT-TSH are tested here to determine the most effective way to reconstruct this quantity and obtain reliable results.

Due to the CT nature of the algorithm, we parallelized our implementation using the Message Passing Interface (MPI) [48]. Despite the additional computation required for the quantum momentum, the runtime of our parallelized implementation is similar to that of a standard TSH approach, as all other calculations are performed in parallel (see discussion at the end of Sec. 5).

In our study, we propose four variants of the implementation of the CT-TSH algorithm, each differing mainly in the method used to evaluate the quantum momentum, and in particular in the approach to determine the appropriate width of the Gaussian functions, as well as in the criteria used to select which atoms are considered in reconstructing the nuclear density.

During the dynamics, the status of each trajectory can change in time as follows:

1. **Active and coupled:** the trajectory is regularly running and is significantly coupled to other trajectories. If $G^{\alpha\beta}$ is the amplitude of the Gaussian representative of trajectory $\beta$, computed at the position of trajectory $\alpha$, the trajectory $\alpha$ is active and coupled if

$$ \frac{\sum_{\beta(\neq\alpha)} G^{\alpha\beta}}{\sum_{\beta} G^{\alpha\beta}} = 1 - \frac{G^{\alpha\alpha}}{\sum_{\beta} G^{\alpha\beta}} > \lambda \qquad (17) $$

where $\lambda$ is a user-defined threshold, which in this case was set to $\lambda = 10^{-2}$. These trajectories obey the CT-TSH equations.

2. **Active and lonely:** the trajectory is regularly running but does not satisfy the inequality (17). In this case, the CT-TSH time evolution trivially reduces to what one would obtain by integrating the TDSE for the single trajectory (as in a standard TSH approach), with the associated overcoherence problem. We apply the overlap-based decoherence correction (ODC) to a trajectory [21], from the time it switches to the lonely condition onwards.

3. **Inactive:** the trajectory is stopped because of technical problems. Inactive trajectories cease to contribute to the CT-TSH dynamics. To elaborate the simulation results, their state probabilities and current states are kept constant from the time each of them turned inactive.

We test two different methods for computing the Gaussian width used in reconstructing the nuclear density. For each Cartesian coordinate $R_{i,\nu}^{\alpha}$, with $i = x, y, z$, a guess for the Gaussian width is determined as the standard deviation over the ensemble of $N_T$ trajectories scaled according to the number of Cartesian coordinates $N_C$:

$$ \sigma_{i,\nu} = \left[ \frac{\sum_{\alpha=1}^{N_T} (R_{i,\nu}^{\alpha})^2}{N_T} - \left( \frac{\sum_{\alpha=1}^{N_T} R_{i,\nu}^{\alpha}}{N_T} \right)^2 \right]^{1/2} \times (N_C)^{1/2}. \qquad (18) $$

The factor $(N_C)^{1/2}$ is introduced in order to reduce the dependence of the $G^{\alpha\beta}/G^{\alpha\alpha}$ ratios on the number of coordinates, so making the method approximately size-consistent. In the first method, the value of $\sigma_{i,\nu}$ is
computed only at the beginning and kept constant throughout the dynamics, meaning that the Gaussians
are frozen. In the second method, the value of $\sigma_{i,\nu}$ is updated at time $t$ if the RHS of Eq. 18 for the
distribution of trajectories at time $t$ exceeds the current value of $\sigma_{i,\nu}$ by a factor of (at least) 1.5 (a user-
defined parameter). It is important to note that $\sigma_{i,\nu}$ is never reduced. While using frozen Gaussians makes
it easier to understand the impact of the variances on the quantum momentum and on the nonadiabatic
dynamics, we also investigated the possibility of allowing $\sigma_{i,\nu}$ to change. This second choice was motivated
by the observation that, with frozen Gaussians, many trajectories became lonely.

An alternative approach is to exclude from the calculation of the quantum momentum the atoms which
are not directly involved in the process considered (“spectator” atoms). Therefore, we introduce a list of
“active” atoms which will contribute in the quantum momentum calculation. Specifically, the active atoms
in our simulations are four: the nitrogen atoms and the two carbon atoms directly bound to them.

Our implementation of the CT-TSH method can be summarized as follows:

1. CT-TSH-FGAC (Coupled-Trajectory Surface Hopping with Frozen Gaussians and All Coordinates): In this approach, we utilize frozen Gaussians and all coordinates contribute to the quantum momentum.

2. CT-TSH-TGAC (Coupled-Trajectory Surface Hopping with Thawed Gaussians and All Coordinates): In this approach, we employ Gaussians whose variance changes over time, and all coordinates contribute to the quantum momentum.

3. CT-TSH-FGLA (Coupled-Trajectory Surface Hopping with Frozen Gaussians and List of Active Atoms): In this approach, we utilize frozen Gaussians and only active atoms contribute to the quantum momentum.

4. CT-TSH-TGLA (Coupled-Trajectory Surface Hopping with Thawed Gaussians and List of Active Atoms): In this approach, we utilize thawed Gaussians, and only active atoms contribute to the quantum momentum.

# 4 Computational details

We test our implementation of CT-TSH in the developer version of MOPAC2002 software [49] by investigating
the photoisomerization dynamics of TAB. Both $n\pi^*$ (which corresponds to the $S_1$ state in the Franck-
Condon region) and $\pi\pi^*$ (which correspond to the $S_2$ state in the Franck-Condon region) excitations were
considered. For comparison purposes, we use three different dynamics methods: TSH, TSH with overlap-
based decoherence correction (TSH-ODC) and CT-TSH. The simulations starting in $S_1$ are performed with
150 trajectories, and those starting in $S_2$ with 100 trajectories. The dynamics with the different approaches
are performed using identical initial conditions. The energies, forces and electronic wavefunctions have been
computed using the FOMO-CI method with a semiempirical AM1 Hamiltonian, replacing the standard set
of AM1 parameters with a set previously reoptimized for azobenzene by some of us [50]. In particular,
a truncated CI within an active space of 13 MOs and 14 electrons is considered, for a total of 94 Slater determinants. All the calculations are performed with the LD algorithm, with an integration time step of 0.1 fs (both for the nuclear and for the electronic degrees of freedom). In the cases where the ODC correction [21] is included, the value of the Gaussian width was 1 a.u. and minimum overlap threshold $S_{min}$ was 0.005. The trajectories were propagated for 1000 fs.

The initial conditions, i.e., geometries and velocities, are determined using the quantum harmonic Wigner sampling technique. To reduce the initial kinetic energy, the C-H stretchings of the phenyl groups were kept frozen during the sampling. Two distinct excitation energy windows are considered: one associated with $n\pi^*$ excitation and the other with $\pi\pi^*$ excitation. For more details, please refer to Table 1. The sampling was performed taking into account the radiative (dipole) transition probability, according to the method outlined in ref. [2].

For the trajectories starting in the $S_1$ state, the $S_1$ classical population is fitted by a delayed exponential function $F_{n\pi^*}$

$$F_{n\pi^*}(t) = \begin{cases} 1 & \text{for } t < t_0 \\ e^{-(t-t_0)/\tau_1} & \text{for } t > t_0 \end{cases} \qquad (19)$$

where $t_0$ is a delay time. The overall lifetime of the $S_1$ state is defined as $\tau_1 = \tilde{\tau}_1 + t_0$. The lifetimes of the $n\pi^*$ ($S_1$) and $\pi\pi^*$ ($S_2$) states for trajectories starting in the $\pi\pi^*$ state are determined using a two-step irreversible kinetic model. To this aim, we fit the $n\pi^*$ population ($F_{n\pi^*}$) and $\pi\pi^*$ population ($F_{\pi\pi^*}$) as

$$F_{\pi\pi^*}(t) = e^{-t/\tau_2} \qquad (20)$$

$$F_{n\pi^*}(t) = \frac{\tau_1}{\tau_1 - \tau_2} (e^{-t/\tau_1} - e^{-t/\tau_2}) \qquad (21)$$

where $\tau_2$ is the $\pi\pi^*$ state lifetime and $\tau_1$ is the $n\pi^*$ state lifetime.

Table 1 summarizes the simulation methods used and provides some information about the results that will be presented in Sec. 5.

The partial photoisomerization quantum yields ($\Phi$) shown in Table 1 were calculated by considering only the trajectories present in the ground state at the end of the simulations ($t = 1000$ fs). Trajectories are labeled as cis if their corresponding CNNC dihedral are less than 45°. To quantify $\Phi$, we define it as the ratio of cis trajectories to the total number of trajectories in the ground state:

$$\Phi = \frac{\text{Number of cis trajectories in the ground state at } t = 1000 \text{ fs}}{\text{Total number of trajectories in the ground state at } t = 1000 \text{ fs}} \qquad (22)$$
<table><caption>Table 1: Nonadiabatic dynamics simulations for TAB: used simulation method, excitation energy windows (eV), fraction of active and lonely trajectories ($F_l$) at the end of the simulated dynamics, fraction of inactive trajectories ($F_{in}$) at the end of the simulated dynamics, fraction of trajectories found in the ground state ($F_{gs}$) at the end of the simulated dynamics, partial photoisomerization quantum yields ($\Phi$), lifetimes $\tau_1$, $\tau_2$ and the delay time $t_0$ in fs.</caption><thead><tr><th>method</th><th>excitation window</th><th>F<sub>l</sub></th><th>F<sub>in</sub></th><th>F<sub>gs</sub></th><th>Φ<sup>a</sup></th><th>τ<sub>1</sub></th><th>τ<sub>2</sub></th><th>t<sub>0</sub></th></tr></thead><tbody><tr><td rowspan="2">CT-TSH-FGAC</td><td>nπ* [2.30 eV, 3.30 eV]</td><td>0.77</td><td>0.01</td><td>0.88</td><td>0.22 ± 0.04</td><td>662</td><td>-</td><td>272</td></tr><tr><td>ππ*[3.80 eV, 4.60 eV]</td><td>0.78</td><td>0.11</td><td>0.81</td><td>0.11 ± 0.04</td><td>463</td><td>99</td><td>-</td></tr><tr><td rowspan="2">CT-TSH-FGLA</td><td>nπ* [2.30 eV, 3.30 eV]</td><td>0.45</td><td>0.01</td><td>0.85</td><td>0.19 ± 0.03</td><td>675</td><td>-</td><td>270</td></tr><tr><td>ππ*[3.80 eV, 4.60 eV]</td><td>0.72</td><td>0.10</td><td>0.74</td><td>0.19 ± 0.05</td><td>291</td><td>384</td><td>-</td></tr><tr><td rowspan="2">CT-TSH-TGAC</td><td>nπ* [2.30 eV, 3.30 eV]</td><td>0.01</td><td>0.03</td><td>0.80</td><td>0.23 ± 0.04</td><td>686</td><td>-</td><td>252</td></tr><tr><td>ππ* [3.80 eV, 4.60 eV]</td><td>0.00</td><td>0.12</td><td>0.71</td><td>0.13 ± 0.04</td><td>443</td><td>131</td><td>-</td></tr><tr><td rowspan="2">CT-TSH-TGLA</td><td>nπ* [2.30 eV, 3.30 eV]</td><td>0.01</td><td>0.01</td><td>0.81</td><td>0.20 ± 0.04</td><td>716</td><td>-</td><td>223</td></tr><tr><td>ππ* [3.80 eV, 4.60 eV]</td><td>0.00</td><td>0.23</td><td>0.39</td><td>0.38 ± 0.08</td><td>398</td><td>738</td><td>-</td></tr><tr><td rowspan="2">TSH</td><td>nπ* [2.30 eV, 3.30 eV]</td><td>-</td><td>0.01</td><td>0.73</td><td>0.18 ± 0.04</td><td>759</td><td>-</td><td>250</td></tr><tr><td>ππ*[3.80 eV, 4.60 eV]</td><td>-</td><td>0.20</td><td>0.46</td><td>0.17 ± 0.05</td><td>476</td><td>435</td><td>-</td></tr><tr><td rowspan="2">TSH-ODC</td><td>nπ* [2.30 eV, 3.30 eV]</td><td>-</td><td>0.01</td><td>0.95</td><td>0.22 ± 0.03</td><td>591</td><td>-</td><td>308</td></tr><tr><td>ππ*[3.80 eV, 4.60 eV]</td><td>-</td><td>0.14</td><td>0.71</td><td>0.13 ± 0.04</td><td>581</td><td>165</td><td>-</td></tr></tbody></table>

# 5 Results

It is well-known that the minimum energy path on the $S_1$ PES of azobenzene connects the *trans* (and the *cis*) geometry to the $S_1/S_0$ conical intersection, located at around 95° of torsion of the CNNC dihedral. As a consequence, the nonadiabatic transitions to $S_0$ in isolated TAB occur mostly in this region or slightly before (CNNC between 90 and 120°). During the decay to the ground state, TAB can isomerize to *cis*-azobenzene (CAB). In this respect, the most important coordinates are the CNNC dihedral and the NNC angles, which is the reason for our choice of the “active” atoms in FGLA and TGLA simulations.

Photoisomerization quantum yields for TAB in low viscosity solvents have been measured by various groups [51, 52, 53, 54, 55, 56, 57, 58], yielding values in the range of 0.20-0.32 for the $n\pi^*$ excitation and 0.09-0.16 for the $\pi\pi^*$ excitation. We present in Table 1, the *trans→cis* partial photoisomerization quantum yields, $\Phi$, for all the methods considered in this work. It is worth noting that trajectories with faster decay tend to predominantly remain in the *trans* side after relaxing to the ground state, while slower trajectories have a higher probability of undergoing photoisomerization. Consequently, we expect a higher photoisomerization quantum yield when all the trajectories reach the ground state. This is why we refer to

<footnote>a Partial quantum yield of the *trans* → *cis* photoisomerization. The binomial standard deviation, obtained as $\sqrt{\Phi(1 - \Phi)/N_T}$, is also shown.</footnote>
the quantum yields presented in Table 1 as “partial”. The calculated Φ values for the trajectories starting
in the $n\pi^*$ state exhibit overall agreement with each other and with the experimental values. We remind
that, besides the simulation method, other factors affect the results and particularly the quantum yields:
possible inaccuracies of the PESs [50], the interactions with the solvent [42, 26], and the sampling of initial
conditions [43]. This last factor, together with the presence of undecayed trajectories after 1 ps, seems to
have a major effect, since previous TSH simulations for the isolated molecule, with decoherence corrections
and Boltzmann sampling of the initial conditions, yielded rather different results: by $n\pi^*$ excitation the
quantum yield was 0.33 [42] and by $\pi\pi^*$ excitation it was 0.20 [26].

The lower photoisomerization quantum yield after $\pi\pi^*$ excitation has been previously discussed [59, 60, 61]. It is related to an exception to Kasha's rule, resulting from a competition between “reactive” and “unreactive” internal conversion (IC), i.e., between IC with photoisomerization or without photoisomerization, respectively. Both processes require a certain degree of progress along the CNNC torsional coordinate. However, the unreactive $S_1 \rightarrow S_0$ IC process can occur earlier (i.e., farther from the 90° midpoint of the torsional pathway) if more vibrational energy is available, as is typically the case with $\pi\pi^*$ excitation. In this context, TSH-ODC, CT-TSH-FGAC, and CT-TSH-TGAC successfully captured this feature, whereas TSH, CT-TSH-FGLA, and CT-TSH-TGLA were unable to do so. The last three approaches showed higher or equivalent photoisomerization quantum yields for $\pi\pi^*$ and $n\pi^*$ excitation.

Fig. 2 and 3 show the state populations as function of time of the trajectories starting in the $n\pi^*$ state,
obtained with the four flavors of CT-TSH (Fig. 2), with TSH (Fig. 3 left) and TSH-ODC (Fig. 3 right).
We remind that the classical populations $F_m(t)$ (thick lines) are calculated as $F_m(t) = N_m(t)/N_T$, where
$N_m(t)$ is the number of trajectories for which the active state at time $t$ is $m$, while the electronic populations
$P_m(t)$ are the averages over all trajectories of the state probabilities $|C_m^\alpha(t)|^2$. In general $F_m(t) \neq P_m(t)$,
and the difference is particularly relevant if decoherence effects are important and not accounted for, like in
TSH, as shown on the left panel of Fig. 3. We note a remarkably good agreement across implementations of
CT-TSH with respect to the electronic populations. However, the thawed Gaussian implementations seems
to lose accuracy at long times if one considers the internal consistency. The partial failure of the thawed
Gaussian option is due to the excessive increase of the $\sigma_{i,\nu}$ widths in time. This is shown in Fig. 4, where
we plot the mass-weighted average $\Sigma = \frac{\sum_{i=1}^{N_C} \sum_{\nu=1}^{N_n} M_{i,\nu}^{1/4} \sigma_{i,\nu}}{N_C}$. Increasing the $\sigma_{i,\nu}$ widths to account for the
broadening of the trajectories distribution as the dynamics proceeds makes the Gaussians very flat: as a
consequence the quantum momentum becomes smaller, because the gradient of the nuclear density almost
vanishes. In such conditions, the CT-TSH method tends to TSH. The complementary drawback of the frozen
Gaussian option, i.e., the large fraction of lonely trajectories, is effectively fixed by the application of ODC
to such trajectories. Note that, before the fraction of lonely trajectories starts to be relevant, i.e., up to
500-600 fs (see Fig. 9), the CT-TSH-FGAC and CT-TSH-FGLA algorithms are able to maintain a good
agreement between $P_m(t)$ and $F_m(t)$. The same is true for the $\pi\pi^*$ excitation, up to 200-300 fs (see Fig. 10).
These considerations confirm the crucial role of the quantum momentum contribution in preserving internal
<figure><img src="image_3.png" /><figcaption>Figure 2: Time evolution of the state populations after $n\pi^*$ excitation for the CT-TSH trajectories. The thick lines represent the classical population, which is the fraction of trajectories running on states $S_0$ (red) and $S_1$ (blue), while the thin lines represent the quantum population, which is the average probability of the $S_0$ and $S_1$ states. The panels report the results of the four implementations of CT-TSH discussed in Sec. 3, as indicated by the acronyms (CT-TSH-FGAC, CT-TSH-FGLA, CT-TSH-TGAC and CT-TSH-TGLA).</figcaption></figure>

consistency throughout the dynamics. Decoherence corrections, such as the ODC algorithm, are designed to restore the internal consistency of the TSH procedure in an *ad hoc* manner, resulting in $F_m(t) \approx P_m(t)$, as seen on the right panel of Fig. 3. As shown in Table 1, for the trajectories starting in the $S_1$ ($n\pi^*$) state, the excited state population decay obtained with TSH is the slowest; the decoherence correction (ODC) changes TSH results, perhaps too much, accelerating the decay significantly; finally, the CT-TSH decay is faster than in TSH, but not as fast as TSH-ODC decay. We can conclude that the decoherence correction applied on TSH, and the coupling among the trajectories, correct TSH in the same direction, accelerating the decay rate.

The faster decay of TSH-ODC was investigated by analyzing the number of hops and “back” hops between the ground state and the first excited state. TSH-ODC exhibited significantly fewer hops (i.e., from $S_1$ to
<figure><img src="image_4.png" /><figcaption>Figure 3: Time evolution of the state populations after $n\pi^*$ excitation for the TSH (left) and TSH-ODC (right) trajectories. The thick lines represent the classical populations, while the thin lines represent the quantum populations. The color code is the same as in Fig. 2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: The mass-weighted average width of the Gaussians over time for thawed Gaussians approaches for the trajectories starting in the $n\pi^*$ state. Right panel: CT-TSH-TGAC. Left panel: CT-TSH-TGLA.</figcaption></figure>

$S_0$) and “back” hops (i.e., from $S_0$ to $S_1$) when compared to other methodologies. Notably, even when the energy difference between the ground state and the first excited state was considerably large ($\ge$ 2.0 eV), a substantial number of hops and “back” hops still occurred in the other methodologies, while they were practically absent in TSH-ODC. Fig. 5 shows the state probabilities of the $S_1$ state after $n\pi^*$ excitation, encompassing all methodologies investigated within this study: the state probabilities in TSH-ODC undergo a sudden rescaling following a hop event, rapidly making the state probability for the new active state very close to one. In this way, the fewest switches hopping probability away from the active state is reduced after a hop in TSH-ODC with respect to CT-TSH.

Fig. 6 and 7 show the state populations as functions of time for the trajectories starting in the $\pi\pi^*$ state, obtained with CT-TSH, TSH and TSH-ODC. It is worth emphasizing that the addition of an extra
<figure><img src="image_6.png" /><figcaption>Figure 5: State probabilities of the $S_1$ state for 10 single trajectories after $n\pi^*$ excitation.</figcaption></figure>

state can significantly increase the complexity of the system. However, our results demonstrate a remarkable agreement between the quantum and classical populations in CT-TSH, by contrast to the TSH approach.

By analyzing Fig. 6 and 7 it is possible to see that the CT-TSH approaches, where all the atoms were considered in the computation of the quantum momentum (CT-TSH-FGAC and CT-TSH-TGCA), presented better agreement with TSH-ODC concerning the populations decay and $\Phi$. The much slower decay of the $S_2$ state observed in CT-TSH when only the CNNC atoms were taken into account in the computation of the quantum momentum (CT-TSH-FGLA and CT-TSH-TGLA) can be attributed to the even higher occurrence of hops and “back” hops between the $S_1$ and $S_2$ states when compared to the other methodologies. For instance, in the case of CT-TSH-TGLA, the occurrence of hops and “back” hops between those states increases by 51% and 61%, respectively, when compared to CT-TSH-FGAC, and by 46% and 36% when compared to CT-TSH-TGAC. It is expected that large values of the Gaussian widths lead to an underestimation of the quantum momentum. By contrast, widths that are too small have the effect to produce a very “bumpy” density: this could potentially lead to a large quantum momentum that changes its direction for small nuclear displacements. Consequently, the *ct* contribution to the time derivative of the electronic coefficients undergoes frequent phase changes, inducing fluctuations in the state probabilities. Then, “back” hops will be more frequent with smaller values of the Gaussian widths, ultimately extending the lifetimes. Notably, when only active atoms are considered in the calculation of the quantum momentum, the values of the widths are 3-4 times smaller than corresponding values when all coordinates are used, both in the frozen-Gaussian and in thawed-Gaussian options (see Fig. 4). In particular, upon analysis of Fig. 8,
it becomes evident that the CT-TSH-TGLA approach displays significantly more frequent fluctuations in the state probabilities compared to the other methodologies, which contributes to an increased occurrence of “back” hops, and, as a result, a longer lifetime of the $S_2$ state is observed.

<figure><img src="image_7.png" /><figcaption>Figure 6: Time evolution of the state populations after $\pi\pi^*$ excitation for the CT-TSH trajectories. The thick lines represent the classical populations, while the thin lines represent the quantum populations. The population of $S_0$ is given in red, of $S_1$ in blue and of $S_2$ in orange.</figcaption></figure>

As shown in Fig. 9 and 10, a significant number of trajectories became lonely when frozen Gaussians are used to compute the quantum momentum. In particular, when considering all atoms for the computation of the quantum momentum, approximately 77% and 78% of trajectories became lonely by the end of the dynamics when the trajectories start on the $n\pi^*$ and $\pi\pi^*$ states, respectively. However, when we limit the active atoms to the carbon and nitrogen atoms bridging the phenyl groups we observe a decrease of the fraction of lonely trajectories to 45% for the trajectories started in the $n\pi^*$ state, and to 72% for the trajectories started in the $\pi\pi^*$ state. This suggests that utilizing fewer coordinates could lead to improved stability of the coupled-trajectory dynamics, because it avoids the involvement of “spectator” atoms in the calculation of the $G^{\alpha\beta}$ factors. Moreover, Fig. 9 and 10 show that all trajectories remain mostly coupled
<figure><img src="image_8.png" /><figcaption>Figure 7: Time evolution of the state populations after $\pi\pi^*$ excitation for the TSH (left) and TSH-ODC (right) trajectories. The thick lines represent the classical populations, while the thin lines represent the quantum populations. The color code is the same as in Fig. 6.</figcaption></figure>

until 400 fs, and only after 500 fs, the number of lonely trajectories begins to increase significantly.

Finally, we investigate the computational time required for propagating CT-TSH trajectories versus TSH trajectories. Table 2 provides information on the total time required for propagating 10 up to 30 trajectories in parallel using CT-TSH and TSH algorithms. In TSH, the total time was calculated as the sum of the time used in each trajectory. For CT-TSH and TSH-S, the trajectories run synchronized in parallel using the same MPI strategy. Consequently, the total time is just the time of one trajectory multiplied by the number of trajectories. The computations were carried out on a machine equipped with an Intel Xeon CPU E5-2450 0 with 32 CPUs, 2 sockets, 8 cores per package, and 2 threads per core.

<table><caption>Table 2: Total time required (in seconds) to propagate 10-30 trajectories with CT-TSH and TSH during 50 fs (500 time steps). The trajectories were started on the $S_1$ ($n\pi^*$) state.</caption><thead><tr><th>$N_T$</th><th>CT-TSH</th><th>TSH</th><th>TSH-S<sup>a</sup></th></tr></thead><tbody><tr><td>10</td><td>3870.4</td><td>1705.4</td><td>2224.6</td></tr><tr><td>20</td><td>13640.2</td><td>4239.2</td><td>7842.8</td></tr><tr><td>30</td><td>19815.9</td><td>8899.7</td><td>11452.5</td></tr></tbody></table>

<sup>a</sup> TSH-S refers for the TSH with synchronized trajectories.

In general, the computational time required in the CT-TSH algorithm is approximately two to three times higher than that of the traditional TSH algorithm. There are several factors that contribute to this increased computational time. One crucial factor to consider is the synchronization of trajectories. In the case of CT-TSH, the propagation of electronic coefficients requires the reconstruction of the nuclear density at each time step. Consequently, the trajectories need to be propagated in a synchronized manner. It is important to stress that the timing of a single time step in the different trajectories can be variable, because of iterative
<figure><img src="image_9.png" /><figcaption>Figure 8: State probabilities of the $S_2$ state for 10 single trajectories after $\pi\pi^*$ excitation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Status of the trajectories as a function of time for CT-TSH starting with the $n\pi^*$ excitation. Right panel: CT-TSH-FGAC. Left panel: CT-TSH-FGLA.</figcaption></figure>

processes (such as SCF) that can require an almost arbitrary number of cycles, which may depend on the nuclear geometry. As a result, all other trajectories must also accommodate these additional computational costs, resulting in a cumulative increase in computational time throughout the dynamics. To evaluate the impact of trajectory synchronization on the overall computational time, we conducted simulations using a synchronized TSH approach, referred to as TSH-S in Table 2. The results clearly demonstrate the crucial role of synchronization in the additional time required in CT-TSH when compared with TSH. However, the most significant time-consuming factor in CT-TSH is the communication among processes necessary
<figure><img src="image_11.png" /><figcaption>Figure 10: Status of the trajectories as a function of time for the CT-TSH-FGAC implementation starting with the $\pi\pi^*$ excitation.</figcaption></figure>

for the computation of quantum momentum. It is important to note that the electronic structure method used in this study is extremely efficient, meaning that the time spent on communication and synchronizing trajectories becomes more significant in comparison. If a more expensive electronic structure method, such as any *ab initio* approach, were employed, the higher computational time required to perform CT-TSH simulations would be significantly less pronounced.

# 6 Conclusion and future perspectives

We presented the implementation and application of a coupled-trajectory scheme for molecular-dynamics simulations of excited-state processes derived from the exact factorization. Specifically, we combined the CT-TSH method with the semiempirical FOMO-CI electronic structure method to investigate the photoisomerization dynamics of *trans*-azobenzene.

To evaluate the efficiency of our approach, we compared the results obtained using the CT-TSH method with the TSH method, with and without ODC decoherence correction. The comparisons included simulations of both $n\pi^*$ and $\pi\pi^*$ excitations. In addition, we tested four implementations of the CT-TSH algorithm, each of them with a distinct method for evaluating the nuclear density used in the computation of the quantum momentum. The nuclear density is represented as a sum of Gaussian functions, each Gaussian being associated with one trajectory. The four variants of the CT-TSH method concern the set of nuclear coordinates considered as arguments of the nuclear density and two different ways to determine the widths of the Gaussian functions.

Our findings showed that the CT-TSH method improved the internal consistency compared to the TSH method and accelerated the excited-state population decay. Notably, using frozen Gaussians (fixed widths) during the dynamics leads to a significant improvement in the agreement between classical and quantum
populations. Internal consistency is lost at long times when thawed Gaussians (variable widths) are used to
calculate quantum momentum. In fact, as the trajectories delocalize in space, the widths of the individual
Gaussians used to reconstruct the nuclear density increases. Therefore, the quantum momentum, i.e., the
spatial derivative of the density, becomes smaller, leading to a decrease of the internal consistency. However,
we believe that by significantly increasing the number of coupled trajectories, we can achieve an accurate
representation of the nuclear distribution in configuration space without requiring a substantial increase
in the Gaussians width. In an ideal scenario, the total number of trajectories would be chosen such that
they do not become lonely when using frozen Gaussians, and the Gaussian width can be kept constant
or almost so. This would ensure that both frozen and thawed Gaussians yield identical results. In our
current work, we have utilized approximately 2 trajectories per number of internal coordinates (TAB has 72
internal coordinates, and we employed 100-150 trajectories). Unfortunately, it was not possible to increase
the number of trajectories due to technical limitations.

The difference in computational time between CT-TSH and TSH is attributed to the need for trajectory
synchronization and to the communication among the processes required for the computation of the quantum
momentum. The relative importance of these two factors depends on the efficiency of the electronic structure
method employed.

Overall, our results demonstrate the good performance of the CT-TSH method and its potential for
simulating excited-state processes, making it a promising tool for broader utilization in the future, for
medium-to-large systems. Moreover, these findings suggest that the CT-TSH method has great potential for
integration with other electronic structure methods, such as *ab initio* approaches.

# 7 Acknowledgments

This work was supported by a public grant from the "Laboratoire d'Excellence Physics Atoms Light Mater" (LabEx PALM) overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference: ANR-10-LABX-0039-PALM), by the ANR Q-DeLight project, Grant No. ANR-20-CE29-0014 of the French Agence Nationale de la Recherche, and by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program, through the "ADI 2019" project funded by the IDEX Paris-Saclay (No. ANR-11-IDEX-0003-02).
# A RK4 method for the numerical integration of the electronic coefficients propagations in CT-TSH using the LD algorithm.

## A.1 Local diabatization algorithm

The basic idea of the LD scheme is to resort to a "locally diabatic" representation, i.e., to a set of electronic states which are specifically diabatic along the nuclear trajectory under consideration. By definition, the diabatic basis spans the same subspace of electronic states and is connected with the adiabatic one by a unitary transformation

$$|\psi\rangle = |\eta\rangle \mathbf{T} \qquad (23)$$

If $\mathbf{H}$ is the Hamiltonian in the diabatic basis ($H_{ij} = \langle\eta_i|\hat{\mathbf{H}}_{\mathbf{B}\mathbf{O}}|\eta_j\rangle$) and $\mathbf{E}$ the diagonal matrix of the electronic energies, we have

$$\mathbf{H}\mathbf{T} = \mathbf{T}\mathbf{E}. \qquad (24)$$

The expansion of the time-dependent wavefunction is

$$\Psi_{el}(t) = \sum_{i} C_{m}(i) |\psi_{m}\rangle = \sum_{i} D_{i}(t) |\eta_{i}\rangle \qquad (25)$$

or

$$\Psi_{el}(t) = |\psi\rangle \mathbf{C}(t) = |\eta\rangle \mathbf{D}(t) \qquad (26)$$

with

$$\mathbf{D}(t) = \mathbf{T}(t)\mathbf{C}(t) \qquad (27)$$

where $\mathbf{D}$ and $\mathbf{C}$ are the vectors of the diabatic and adiabatic coefficients, respectively. The diabatic basis, in the LD scheme, is redefined at each time step. At the beginning of the time step ($t=0$), the transformation matrix is chosen to be the identity matrix, thus $\eta(0) = \psi(0)$. By choosing the diabatic states to be as invariant as possible within the time step considered, the dynamic couplings $\langle\eta_i|\frac{\partial}{\partial t}|\eta_j\rangle$ can be neglected. However, the coupling vector vanishes only along the advancement coordinate identified by the velocity vector $\mathbf{v}$: in this sense the $\eta$ are "locally" diabatic for the given trajectory.

If we just consider the TDSE for independent trajectories, as in TSH, we can easily obtain the time-evolution of a set of diabatic coefficients we shall call $\mathbf{D}_{tsh}$. By inserting the electronic wavefunction expansion 26 into the TDSE (Eq. 2) we obtain

$$\dot{\mathbf{D}}_{tsh} = -\frac{i}{\hbar}\mathbf{H}(t)\mathbf{D}_{tsh}. \qquad (28)$$

Within the LD algorithm, the transformation matrix at the end of the time step, $\mathbf{T}(\Delta t)$, is obtained by Löwdin orthonormalization of the wavefunction overlap matrix

$$S_{ml}(\Delta t) = \langle\psi_{m}(0)|\psi_{l}(\Delta t)\rangle. \qquad (29)$$
Hence, from $T(\Delta t)$ one gets $H(\Delta t)$

$$H(\Delta t) = T(\Delta t)E(\Delta t)T^t(\Delta t). \qquad (30)$$

The matrix $H(t)$ at intermediate times in the interval $[0, \Delta t]$ is obtained by linear interpolation

$$H(t) = E(0) + [H(\Delta t) - E(0)] \frac{t}{\Delta t}. \qquad (31)$$

The essential point is that, unlike the NAC vectors, the diabatic quantities (such as $H$) depend smoothly on the nuclear coordinates and can be easily and accurately interpolated. Eq. 28 can be easily solved for any time $t \in [0, \Delta t]$.

To obtain the electronic evolution equation for the CT-TSH method, we can combine Eq. 28 with Eq. 11. This yields an expression equivalent to Eq. 5 in diabatic basis

$$\dot{D}(t) = -\frac{i}{\hbar} H(t) D_{tsh} + \dot{D}_{ct}(t) \qquad (32)$$

where the elements of $\dot{D}_{ct}(t)$ can be written as

$$\dot{D}_{m,ct}(t) = \sum_{\nu=1}^{N_C} \frac{\mathcal{P}_{\nu}}{\hbar M_{\nu}} \cdot \left( \mathbf{f}_{\nu,m} - \sum_{l} |C_l(t)|^2 \mathbf{f}_{\nu,l} \right) D_m(t). \qquad (33)$$

Eq. 33 was obtained by combining Eq. 13 with the vector potential approximation obtained from the accumulated forces, as shown in Eq. 14. It should be noted that $A_{\nu} \approx \sum_{l} |C_l(t)|^2 \mathbf{f}_{\nu,l}$, and $A_{\nu}$ remains invariant under changes of basis. This is the reason why we retain the adiabatic coefficients in our expression. Additionally, we introduce an approximation by assuming that the accumulated adiabatic forces $\mathbf{f}_{\nu,m}$ and the quantum momentum $\mathcal{P}_{\nu}$ remain constant within the integration time step. This simplification enables us to streamline the expression and avoid the diabatization of the accumulated adiabatic forces, since in the LD scheme, at the beginning of the time step, $D(0) \equiv A(0)$.

Finally, the adiabatic coefficients can be recovered by inverting Eq. 27:

$$C(t) = T^t(t)D(t). \qquad (34)$$

## A.2 RK4 algorithm

To proceed, we will evaluate equation 33 by employing the RK4 algorithm. In this particular case, the function to be solved corresponds to $\dot{D}(t, D)$. It is worth noting that at the beginning of each time step, we initialize the diabatic coefficients to be equal to the adiabatic coefficients, denoted as $C(0) = D(0)$. Additionally, we interpolate the diabatic Hamiltonian using equation 31. With these initial conditions established, we can determine the slope ($k_1, k_2, k_3$ and $k_4$) through the following procedure (from this point onward, all expressions are expressed in atomic units):
### A.2.1 $k_1$ in LD

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

$$H(0) = E(0) \qquad (35)$$

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

$$D(0) = C(0) \qquad (36)$$

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

$$\dot{D}_{m,ct}(0) = \sum_{\nu=1}^{N_C} \frac{\mathcal{P}_{\nu}}{M_{\nu}} \cdot \left( \mathbf{f}_{\nu,m} - \sum_{l} |C_l(0)|^2 \mathbf{f}_{\nu,l} \right) D_m(0) \qquad (37)$$

$$k_1 = -i\mathbf{E}(0)\mathbf{C}(0) + \dot{D}_{ct}(0) \qquad (38)$$

### A.2.2 $k_2$ in LD

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

$$H(\Delta t/2) = E(0) + \frac{1}{2}[H(\Delta t) - E(0)] \qquad (39)$$

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

$$D(\Delta t/2) = C(0) + \frac{(\Delta t \cdot k_1)}{2} \qquad (40)$$

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

$$\dot{D}_{m,ct}(\Delta t/2) = \sum_{\nu=1}^{N_C} \frac{\mathcal{P}_{\nu}}{M_{\nu}} \cdot \left( \mathbf{f}_{\nu,m} - \sum_{l} \left| \sum_{j} T_{jl}(\Delta t/2) D_j(\Delta t/2) \right|^2 \mathbf{f}_{\nu,l} \right) D_j(\Delta t/2) \qquad (41)$$

$$k_2 = -i \left[ \mathbf{E}(0) + \frac{1}{2}[\mathbf{H}(\Delta t) - \mathbf{E}(0)] \right] \left[ \mathbf{C}(0) + \frac{(\Delta t \cdot k_1)}{2} \right] + \dot{D}_{ct}(\Delta t/2) \qquad (42)$$

We know that at $t = 0$, $T(0) = \mathbf{1}$ and at $t = \Delta t$, $T(\Delta t) = \mathbf{S}(\Delta t)$ (as shown in Eq. 29). The matrix $T(\Delta t/2)$ is estimated by linear interpolation. A proper way to obtain $T(\Delta t/2)$ would be the diagonalization of $H(\Delta t/2)$. However, this introduces a sign arbitrariness which would require careful resolution (not implemented in this first setup of the CT-TSH method).

### A.2.3 $k_3$ in LD

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

$$H(\Delta t/2) = E(0) + \frac{1}{2}[H(\Delta t) - E(0)] \qquad (43)$$

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

$$D(\Delta t/2) = C(0) + \frac{(\Delta t \cdot k_2)}{2} \qquad (44)$$

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

$$\dot{D}_{m,ct}(\Delta t/2) = \sum_{\nu=1}^{N_C} \frac{\mathcal{P}_{\nu}}{\hbar M_{\nu}} \left( \mathbf{f}_{\nu,m} - \sum_{l} \left| \sum_{j} T_{jl}(\Delta t/2) D_j(\Delta t/2) \right|^2 \mathbf{f}_{\nu,l} \right) D_j(\Delta t/2) \qquad (45)$$
$$k_3 = -i \left[ \mathbf{E}(0) + \frac{1}{2} [\mathbf{H}(\Delta t) - \mathbf{E}(0)] \right] \left[ \mathbf{C}(0) + \frac{(\Delta t \cdot k_2)}{2} \right] + \dot{\mathbf{D}}_{ct}(\Delta t/2) \quad (46)$$

### A.3 $k_4$ in LD

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

$$\mathbf{H}(\Delta t) = \mathbf{T}(\Delta t)\mathbf{E}(\Delta t)\mathbf{T}(\Delta t)^t \quad (47)$$

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

$$\mathbf{D}(\Delta t) = \mathbf{C}(0) + \Delta t \cdot k_3 \quad (48)$$

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

$$\dot{\mathbf{D}}_{m,ct}(\Delta t) = \sum_{\nu=1}^{N_C} \frac{\mathcal{P}_{\nu}}{\hbar M_{\nu}} \left( \mathbf{f}_{\nu,m} - \sum_{l} \left| \sum_{j} T_{jl}(\Delta t) D_j(\Delta t) \right|^2 \mathbf{f}_{\nu,l} \right) D_j(\Delta t) \quad (49)$$

$$k_4 = -i \mathbf{H}(\Delta t) [\mathbf{A}(0) + \Delta t \cdot k_3] + \dot{\mathbf{D}}_{ct}(\Delta t) \quad (50)$$

The new adiabatic coefficients are:

$$\mathbf{C}(\Delta t) = \mathbf{T}(\Delta t)^t \left[ \mathbf{C}(0) + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4) \right]. \quad (51)$$

## References

[1] John C. Tully. Mixed quantum-classical dynamics. Faraday Discuss., 110:407–419, 1998.

[2] Maurizio Persico and Giovanni Granucci. An overview of nonadiabatic dynamics simulations methods, with focus on the direct approach versus the fitting of potential energy surfaces. Theor. Chem. Acc., 133(9):1526, 2014.

[3] Elizabeth Brunk and Ursula Rothlisberger. Mixed quantum mechanical/molecular mechanical molecular dynamics simulations of biological systems in ground and electronically excited states. Chem. Rev., 115(12):6217–6263, 2015.

[4] Rachel Crespo-Otero and Mario Barbatti. Recent advances and perspectives on nonadiabatic mixed quantum-classical dynamics. Chem. Rev., 118(15):7026–7068, 2018.

[5] Basile F. E. Curchod and Todd J. Martínez. Ab initio nonadiabatic quantum molecular dynamics. Chem. Rev., 118(7):3305–3336, 2018.

[6] Leticia González and Roland Lindh, editors. Quantum Chemistry and Dynamics of Excited States. Wiley, November 2020.

[7] Maurizio Persico, Giovanni Granucci, and Davide Accomasso. The quantum decoherence problem in nonadiabatic trajectory methods. In Comprehensive Computational Chemistry, Russell Boyd and Manuel Yañez. Elsevier, 2023.
[8] Alexey V. Akimov. Fundamentals of trajectory-based methods for nonadiabatic dynamics. In *Comprehensive Computational Chemistry, Russell Boyd and Manuel Yañez*. Elsevier, 2023.

[9] Eric R. Bittner and Peter J. Rossky. Quantum decoherence in mixed quantum-classical systems: Nonadiabatic processes. *J. Chem. Phys.*, 103(18):8130–8143, 11 1995.

[10] Basile F. E. Curchod and Ivano Tavernelli. On trajectory-based nonadiabatic dynamics: Bohmian dynamics versus trajectory surface hopping. *J. Chem. Phys.*, 138(18):184112, 05 2013.

[11] Basile F. E. Curchod, Federica Agostini, and E. K. U. Gross. An exact factorization perspective on quantum interferences in nonadiabatic dynamics. *J. Chem. Phys.*, 145(3):034103, 07 2016.

[12] Yinan Shu and Donald G. Truhlar. Decoherence and its role in electronically nonadiabatic dynamics. *J. Chem. Theory Comput.*, 19(2):380–395, 2023.

[13] John C. Tully. Molecular dynamics with electronic transitions. *J. Chem. Phys.*, 93(2):1061–1071, 1990.

[14] Ali Abedi, Neepa T. Maitra, and E. K. U. Gross. Exact factorization of the time-dependent electron-nuclear wave function. *Phys. Rev. Lett.*, 105:123002, Sep 2010.

[15] Seung Kyu Min, Federica Agostini, and E. K. U. Gross. Coupled-trajectory quantum-classical approach to electronic decoherence in nonadiabatic processes. *Phys. Rev. Lett.*, 115:073001, Aug 2015.

[16] Lea M. Ibele, Carlotta Pieroni, Francesco Talotta, Basile F.E. Curchod, David Lauvergnat, and Federica Agostini. Exact factorization of the electron-nuclear wavefunction: Fundamentals and algorithms. In *Comprehensive Computational Chemistry, Russell Boyd and Manuel Yañez*. Elsevier, 2023.

[17] Carlotta Pieroni and Federica Agostini. Nonadiabatic dynamics with coupled trajectories. *J. Chem. Theory Comput.*, 17(10):5969–5991, 2021.

[18] Jian-Yun Fang and Sharon Hammes-Schiffer. Improvement of the internal consistency in trajectory surface hopping. *J. Phys. Chem. A*, 103(47):9399–9407, 1999.

[19] Chaoyuan Zhu, Shikha Nangia, Ahren W. Jasper, and Donald G. Truhlar. Coherent switching with decay of mixing: An improved treatment of electronic coherence for non-Born–Oppenheimer trajectories. *J. Chem. Phys.*, 121(16):7658–7670, 2004.

[20] Giovanni Granucci and Maurizio Persico. Critical appraisal of the fewest switches algorithm for surface hopping. *J. Chem. Phys.*, 126(13):134114, 2007.

[21] Giovanni Granucci, Maurizio Persico, and Alberto Zoccante. Including quantum decoherence in surface hopping. *J. Chem. Phys.*, 133(13):134111, 2010.
[22] Joseph E. Subotnik, Wenjun Ouyang, and Brian R. Landry. Can we derive Tully's surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence. *J. Chem. Phys.*, 139(21):214107, 2013.

[23] Joseph E. Subotnik, Amber Jain, Brian Landry, Andrew Petit, Wenjun Ouyang, and Nicole Bellonzi. Understanding the surface hopping view of electronic transitions and decoherence. *Annu. Rev. Phys. Chem.*, 67(1):387-417, 2016.

[24] Patricia Vindel-Zandbergen, Lea M. Ibele, Jong-Kwon Ha, Seung Kyu Min, Basile F. E. Curchod, and Neepa T. Maitra. Study of the decoherence correction derived from the exact factorization approach for nonadiabatic dynamics. *J. Chem. Theory Comput.*, 17(7):3852-3862, 2021.

[25] Giovanni Granucci, Maurizio Persico, and Alessandro Toniolo. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. *J. Chem. Phys.*, 114(24):10608-10615, 2001.

[26] Valentina Cantatore, Giovanni Granucci, and Maurizio Persico. Simulation of the $\pi \rightarrow \pi^*$ photodynamics of azobenzene: Decoherence and solvent effects. *Comput. Theor. Chem.*, 1040-1041:126-135, 2014.

[27] Davide Accomasso, Giovanni Granucci, Meilani Wibowo, and Maurizio Persico. Delocalization effects in singlet fission: Comparing models with two and three interacting molecules. *J. Chem. Phys.*, 152(24):244125, 06 2020.

[28] Eduarda Sangiogo Gil, Maurizio Persico, and Giovanni Granucci. Frenkel exciton photodynamics of self-assembled monolayers of azobiphenyls. *J. Chem. Phys.*, 157(16):161101, 10 2022.

[29] Giacomo Salvadori, Veronica Macaluso, Giulia Pellicci, Lorenzo Cupellini, Giovanni Granucci, and Benedetta Mennucci. Protein control of photochemistry and transient intermediates in phytochromes. *Nat. Commun.*, 13(1), 2022.

[30] Ali Abedi, Neepa T. Maitra, and E. K. U. Gross. Correlated electron-nuclear dynamics: Exact factorization of the molecular wavefunction. *J. Chem. Phys.*, 137(22):22A530, 2012.

[31] Federica Agostini and E. K. Gross. Ultrafast dynamics with the exact factorization. *Eur. Phys. J. B*, 94(9), 2021.

[32] Federica Agostini, Seung Kyu Min, Ali Abedi, and E. K. U. Gross. Quantum-classical nonadiabatic dynamics: Coupled- vs independent-trajectory methods. *J. Chem. Theory Comput.*, 12(5):2127-2143, 2016.

[33] Seung Kyu Min, Federica Agostini, Ivano Tavernelli, and E. K. U. Gross. Ab initio nonadiabatic dynamics with coupled trajectories: A rigorous approach to quantum (de)coherence. *J. Phys. Chem. Lett.*, 8(13):3048-3055, 2017.
[34] Graeme H. Gossel, Federica Agostini, and Neepa T. Maitra. Coupled-trajectory mixed quantum-classical algorithm: A deconstruction. *J. Chem. Theory Comput.*, 14(9):4513-4529, 2018.

[35] In Seong Lee, Jong-Kwon Ha, Daeho Han, Tae In Kim, Sung Wook Moon, and Seung Kyu Min. Pyunixmd: A python-based excited state molecular dynamics package. *J. Comput. Chem.*, 42(24):1755-1766, 2021.

[36] Evaristo Villaseco Arribas, Federica Agostini, and Neepa T. Maitra. Exact factorization adventures: A promising approach for non-bound states. *Molecules*, 27(13), 2022.

[37] Evaristo Villaseco Arribas and Neepa T. Maitra. Energy-conserving coupled trajectory mixed quantum-classical dynamics. *J. Chem. Phys.*, 158(16):161105, 2023.

[38] Jong-Kwon Ha, In Seong Lee, and Seung Kyu Min. Surface hopping dynamics beyond nonadiabatic couplings for quantum coherence. *J. Phys. Chem. Lett.*, 9(5):1097-1104, 2018.

[39] Jong-Kwon Ha and Seung Kyu Min. Independent trajectory mixed quantum-classical approaches based on the exact factorization. *J. Chem. Phys.*, 156(17):174109, 2022.

[40] Daeho Han, Jong-Kwon Ha, and Seung Kyu Min. Real-space and real-time propagation for correlated electron-nuclear dynamics based on exact factorization. *J. Chem. Theory Comput.*, 19(8):2186-2197, 2023.

[41] Federica Agostini. An exact-factorization perspective on quantum-classical approaches to excited-state dynamics. *Eur. Phys. J. B*, 91(7), 2018.

[42] Teresa Cusati, Giovanni Granucci, and Maurizio Persico. Photodynamics and time-resolved fluorescence of azobenzene in solution: A mixed quantum-classical simulation. *J. Am. Chem. Soc.*, 133(13):5109-5123, 2011.

[43] Carlotta Pieroni, Filippo Becuzzi, Luigi Creatini, Giovanni Granucci, and Maurizio Persico. Effect of initial conditions sampling on surface hopping simulations in the ultrashort and picosecond time range. azomethane photodissociation as a case study. *J. Chem. Theory Comput.*, 19(9):2430-2445, 2023.

[44] Felix Plasser, Giovanni Granucci, Jiri Pittner, Mario Barbatti, Maurizio Persico, and Hans Lischka. Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer. *J. Chem. Phys.*, 137(22):22A514, 2012.

[45] Isabella C. D. Merritt, Denis Jacquemin, and Morgane Vacher. Nonadiabatic coupling in trajectory surface hopping: How approximations impact excited-state reaction dynamics. *J. Chem. Theory Comput.*, 19(6):1827-1842, 2023.
[46] Tian Qiu, Clàudia Climent, and Joseph E. Subtonik. A practical approach to wave function propagation, hopping probabilities, and time steps in surface hopping calculations. *J. Chem. Theory Comput.*, 0(0):null, 2023.

[47] Neus Aguilera-Porta, Inés Corral, Jordi Munoz-Muriedas, and Giovanni Granucci. Excited state dynamics of some nonsteroidal anti-inflammatory drugs: A surface-hopping investigation. *Comput. Theor. Chem.*, 1152:20-27, 2019.

[48] CORPORATE The MPI Forum. Mpi: A message passing interface. In *Proceedings of the 1993 ACM/IEEE Conference on Supercomputing*, Supercomputing '93, page 878-883, New York, NY, USA, 1993. Association for Computing Machinery.

[49] James J. P. Stewart. *MOPAC2002; Fujitsu Limited: Tokyo, Japan*, 2002.

[50] Teresa Cusati, Giovanni Granucci, Emilio Martínez-Núñez, Francesca Martini, Maurizio Persico, and Saulo Vázquez. Semiempirical hamiltonian for simulation of azobenzene photochemistry. *J. Phys. Chem. A*, 116:98, 2012.

[51] Jacques Ronayette, René Arnaud, Philippe Lebourgeois, and Jacques Lemaire. Isomérisation photochimique de l'azobenzène en solution. i. *Can. J. Chem.*, 52(10):1848-1857, 1974.

[52] Pietro Bortolus and Sandra Monti. Cis trans photoisomerization of azobenzene. solvent and triplet donors effects. *J. Phys. Chem.*, 83(6):648-652, 1979.

[53] Günter Gauglitz and Stephan Hubig. Chemical actinometry in the uv by azobenzene in concentrated solution: A convenient method. *J. Photochem.*, 30(2):121-125, 1985.

[54] Narcisse Siampiringue, Ghislain Guyot, Sandra Monti, and Pietro Bortolus. The cis-trans photoisomerization of azobenzene: an experimental re-examination. *J. Photochem.*, 37(1):185-188, 1987.

[55] Arthur W. Adamson, Arnd Vogler, Horst Kunkely, and Rudolf Wachter. Photocalorimetry. enthalpies of photolysis of trans-azobenzene, ferrioxalate and cobaltioxalate ions, chromium hexacarbonyl, and dirhenium decarbonyl. *J. Am. Chem. Soc.*, 100(4):1298-1300, 1978.

[56] Hermann Rau and Erik Lueddecke. On the rotation-inversion controversy on photoisomerization of azobenzenes. experimental proof of inversion. *J. Am. Chem. Soc.*, 104(6):1616-1620, 1982.

[57] Hermann Rau. Further evidence for rotation in the $\pi\pi^*$ and inversion in the $n\pi^*$ photoisomerization of azobenzenes. *J. Photochem.*, 26(2):221-225, 1984.

[58] Vít Ladányi, Pavel Dvořák, Jamaludin Al Anshori, Ľubica Vetráková, Jakob Wirz, and Dominik Heger. Azobenzene photoisomerization quantum yields in methanol redetermined. *Photochem. Photobiol. Sci.*, 16:1757-1761, 2017.
[59] Cosimo Ciminelli, Giovanni Granucci, and Maurizio Persico. The photoisomerization mechanism of azobenzene: A semiclassical simulation of nonadiabatic dynamics. *Chem. Eur. J.*, 10(9):2327–2341, 2004.

[60] Artur Nenov, Rocio Borrego-Varillas, Aurelio Oriana, Lucia Ganzer, Francesco Segatta, Irene Conti, Javier Segarra-Marti, Junko Omachi, Maurizio Dapor, Simone Taioli, Cristian Manzoni, Shaul Mukamel, Giulio Cerullo, and Marco Garavelli. Uv-light-induced vibrational coherences: The key to understand kasha rule violation in trans-azobenzene. *J. Chem. Phys. Lett.*, 9(7):1534–1541, 2018.

[61] Jimmy K. Yu, Christoph Bannwarth, Ruibin Liang, Edward G. Hohenstein, and Todd J. Martínez. Nonadiabatic dynamics simulation of the wavelength-dependent photochemistry of azobenzene excited to the $n\pi^*$ and $\pi\pi^*$ excited states. *J. Am. Chem. Soc.*, 142(49):20680–20690, 2020.