# Velocity Rescaling in Surface Hopping Based on Atomic Contributions to Electronic Transitions

Eduarda Sangiogo-Gil, Lea Ibele, Richard Bleyer, Leticia González

► To cite this version:

Eduarda Sangiogo-Gil, Lea Ibele, Richard Bleyer, Leticia González. Velocity Rescaling in Surface Hopping Based on Atomic Contributions to Electronic Transitions. Journal of Chemical Theory and Computation, 2025, 21 (17), pp.8278-8290. 10.1021/acs.jctc.5c00737. hal-05288243

HAL Id: hal-05288243

https://hal.science/hal-05288243v1

Submitted on 1 Oct 2025

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>

# Velocity Rescaling in Surface Hopping Based on Atomic Contributions to Electronic Transitions

Eduarda Sangiogo-Gil,* Lea M. Ibele, Richard Bleyer, and Leticia González*

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

Cite This: J. Chem. Theory Comput. 2025, 21, 8278–8290

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

Read Online

## ACCESS

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

Metrics & More

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

Article Recommendations

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

Supporting Information

**ABSTRACT:** Surface hopping is a widely used method for simulating nonadiabatic dynamics, in which nuclear motion follows classical trajectories and electronic transitions occur stochastically. To ensure energy conservation during these transitions, atomic velocities must be adjusted. Traditional velocity rescaling methods either apply a uniform adjustment to atomic velocities, which can lead to size-consistency issues, or rely on nonadiabatic coupling vectors, which are computationally expensive and may not always be available. Here, we introduce two novel velocity rescaling methods that incorporate atomic contributions to electronic transitions, derived from the one-electron transition density matrix or the density difference between states for a given transition. The first method, *excitation-weighted velocity rescaling*, redistributes

kinetic energy among atoms proportionally to their contributions to the electronic transition. This is achieved through a weighted scaling factor, computed from the population analysis of the one-electron transition density matrix or the density difference of the two states involved in the transition. The second method, *excitation-thresholded velocity rescaling*, adjusts the velocities only of atoms whose contributions exceed a predefined threshold, preventing unnecessary energy redistribution to atoms with minimal involvement in the excitation. We validate these approaches through excited-state dynamics simulations of fulvene and 1H-1,2,3-triazole. Our results show that excitation-weighted velocity rescaling closely reproduces the adjustments based on nonadiabatic coupling vectors for both fulvene and 1H-1,2,3-triazole.

Atomic contribution to the electronic transition

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

Velocity adjustment

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

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

## 1. INTRODUCTION

The surface hopping (SH) method is a widely used approach for simulating nonadiabatic dynamics, where electrons are treated quantum mechanically and nuclei follow classical trajectories, bridging quantum and classical descriptions of molecular motion.¹⁻⁸ The most commonly used version of SH—based on the fewest-switches algorithm—was introduced by Tully in 1990.⁹ This algorithm operates on the premise that nuclear motion is predominantly adiabatic, with nonadiabatic transitions occurring only momentarily and within localized regions of the configuration space. Tully proposed modeling these transitions as instantaneous “hops” between adiabatic potential energy surfaces (PESs). Within this framework, multiple independent trajectories are simulated, forming a statistical ensemble. The relative fraction of trajectories on each PES is used to approximate the population of each quantum state in a realistic dynamical process. The hops are not arbitrary: their probability is dictated by the electronic coupling between states and becomes significant in regions where the energy difference between PESs is small.¹⁰,¹¹

The widespread success of SH over the past three decades stems from its simplicity and intuitive nature. The concept of nuclei evolving classically while stochastically transitioning between states makes the results easy to interpret. Additionally,

SH is straightforward to implement computationally and, despite certain approximations, often yields reliable results across a variety of systems. One key advantage that makes it very attractive is its efficiency in handling large polyatomic systems, as calculations can be performed on-the-fly, eliminating the need for precomputed PESs and couplings.¹² SH can also benefit from predefined PESs, such as linear vibronic coupling models (LVC), making the propagation in full dimensionality extremely efficient.¹³⁻¹⁵ The latter avoids the nontrivial task of selecting the most relevant modes in wavepacket propagation methods, while preserving high efficiency due to the predefined nature of the PESs. However, LVC models also impose limitations, as the linear coupling assumption around a reference geometry restricts their ability to capture anharmonic or large-amplitude nuclear motions. As a result, SH simulations based on LVCs are more suitable for

Received: May 7, 2025

Revised: August 7, 2025

Accepted: August 18, 2025

Published: August 22, 2025

<figure><img src="image_11.png" /></figure>
studying photophysics (e.g., population dynamics) than photochemistry involving significant bond-breaking or structural rearrangements.

Nonetheless, the advantages of SH come with several approximations and adjustable parameters that can influence the resulting dynamics. Notable examples include the treatment of decoherence¹⁶⁻²⁰ the algorithms for calculating hop probabilities,⁹,²¹,²² zero-point energy leakage correction,²³⁻²⁶ the choice of initial conditions,²⁶,²⁷ and the method for rescaling velocities after a hop.²⁸⁻³³ The latter is the focus of this study. When the system hops to a new state, the kinetic energy of the nuclei must be adjusted to ensure classical energy conservation. This adjustment is necessary because the new state has a different potential energy than the previous one, and the total energy (i.e., kinetic plus potential energy) needs to be constant. Accordingly, the kinetic energy is rescaled by modifying the velocity vector v such that the new kinetic plus potential energy equals the total energy of the previous time step, before the hop. This means that when a transition occurs between electronic states, in SH the velocities must be adjusted to conserve total energy. An approach commonly used is to rescale the velocity along its direction by a factor determined as

$$ \mathbf{v}^{\text{new}} = \mathbf{v}^{\text{old}} \sqrt{\frac{E_{\text{kin}}^{\text{new}}}{E_{\text{kin}}^{\text{old}}} \quad (1)} $$

where $E_{\text{kin}}^{\text{new}} = E_{\text{tot}} - E_{\text{pot}}^{\text{new}}$, with $E_{\text{tot}}$ being the total energy, $E_{\text{pot}}^{\text{new}}$ the potential energy of the new electronic state, and $E_{\text{kin}}^{\text{old}}$ is given by $E_{\text{kin}}^{\text{old}} = E_{\text{tot}} - E_{\text{pot}}^{\text{old}}$, where $E_{\text{pot}}^{\text{old}}$ is the potential energy before the hop.²,⁶ This approach is simple and computationally cheap, as it does not require additional quantities, such as nonadiabatic coupling vectors (NACVs), which are cumbersome to obtain. However, and critically, it introduces a size-consistency problem in the SH algorithm.⁶,³³ The concept of size-consistency is reminiscent of electronic structure theory but can be used in an analogous way for nonadiabatic dynamics. Consider the following two scenarios: (i) a molecule A; and (ii) the same molecule A, but with a molecule B at an infinite distance. In the second case, since A and B do not interact, the dynamics of A should remain identical as in the first case. However, if a hop in A requires more energy than is available in its total energy, the results will differ. In case (i), where only A is present, the hop is “frustrated” and cannot occur. In case (ii), where A and B form a larger system, the total kinetic energy of the system may be sufficient to allow the hop. Here, rescaling the velocity of the entire system is unphysical, as the energy needed for the transition in A is effectively “borrowed” by slowing down atoms in B, despite A and B being noninteracting.

This size-consistency problem is not limited to infinitely separated systems. It also arises in finitely separated systems that contain weakly interacting subsystems,³⁴⁻³⁶ and it is particularly pronounced in large systems, such as chromophores in solution or in biological environments, where not all parts of the system are necessarily involved in the transition. For example, consider a system composed of multiple chromophores (such as those in refs 35 and 37⁻⁴⁰). In such cases, the excitation can be localized in a subset of the chromophores, a single chromophore, or even a specific region of a chromophore. When SH dynamics are performed for these systems using the velocity rescaling presented in eq 1, the total

kinetic energy of all the chromophores is available. This means that even if the excitation is strictly localized to a single chromophore or region, the system as a whole may have sufficient kinetic energy for a transition to occur. This could allow a hop to occur even when it is physically unmotivated, creating an artificial situation in which transitions occur that would not be plausible if the excitation were strictly confined. In simulations of systems with multiple chromophores, such as the 12-azobenzene chromophore setup in ref 35, the system's total kinetic energy may always appear sufficient to allow a hop, even though the excitation is confined to one chromophore.

For simulations of chromophores in solution performed using a hybrid quantum mechanical/molecular mechanical (QM/MM) approach, it is common to consider only the kinetic energy of the atoms in the chromophore, which is typically the one molecule included in the QM region.⁴¹ However, a large QM region can also encompass atoms that do not directly contribute to the electronic transition of interest.

A more systematic and size-consistent method, originally proposed by Tully, involves adjusting the velocity component along the NACV (see Appendix A).²⁸,⁴² He demonstrated the physical motivation for this choice by deriving the Pechukas force—a semiclassical expression for the effective force that governs classical-like trajectory methods,⁴³ such as Ehrenfest dynamics and SH, in a two-state quantum system that includes nonadiabatic effects. This force consists of three components: two aligned along the potential energy gradients of the current and target electronic states, and a third along the NACV. The latter, termed the transition force, is the only component proportional to the potential energy gap, ensuring energy conservation during the nonadiabatic transition.

While it is widely accepted that nuclear momentum (or, equivalently, nuclear velocity, since in surface hopping, momentum and velocity are related by momentum equal mass times the velocity. Since the mass is constant for each particle, rescaling the velocity is equivalent to rescaling the momentum) should be corrected along the NACV direction⁹,²⁹,³³,⁴⁴⁻⁴⁷ the direct computation of NACVs is often impractical. NACVs require the evaluation of wave function derivatives with respect to nuclear coordinates, which can be computationally expensive or even unavailable—as is the case in many electronic structure codes. It is worth noting that, although NACV are formally first-order derivatives of the electronic wave functions, their computation is not as widely available or as straightforward as that of energy gradients. NACVs involve cross-terms between different electronic states, $\langle\psi_l|\nabla_R\psi_j\rangle$, and require careful phase tracking and overlap evaluations. While gradients are routinely available in many electronic structure methods, analytic NACVs are only implemented in a limited set of methods and are often more computationally demanding.

To circumvent this limitation, several SH algorithms avoid explicit NACV calculations when solving the electronic time-dependent Schrödinger equation.⁴⁸⁻⁵⁰ In this context, many SH implementations estimate the time-derivative couplings, $\langle\psi_l|\frac{\partial\psi_j}{\partial t}\rangle = \frac{\partial R}{\partial t}\langle\psi_l|\nabla_R\psi_j\rangle = \mathbf{v} \cdot \mathbf{d}_{ij}$ using numerical differentiation.⁵¹,⁵² A more robust alternative involves the evaluation of wave function overlaps between consecutive time steps. Among these strategies, one of the most widely employed is the local diabatization scheme,⁵³,⁵⁴ which circumvents the
explicit calculation of nonadiabatic coupling vectors by constructing a diabatic-like basis from the overlap matrix of adiabatic states at two successive nuclear geometries. In such cases, the most common procedure for rescaling nuclear velocities after a surface hop follows eq 1.

Recently, Toldo et al.³³ systematically investigated various velocity adjustment schemes in SH dynamics. They assessed the gradient difference method as an approximation for the NACV and explored an alternative approach that reduces the available kinetic energy by distributing it across the system's vibrational degrees of freedom before a nonadiabatic transition. The latter approach rescales velocities along the momentum direction; however, it mitigates the occurrence of so-called "unavoidable" back-hoppings by reducing the kinetic energy available for hops. However, it introduces size-consistency issues, as the uniform distribution of kinetic energy among all vibrational modes can slow down the motion of the degrees of freedom that are not directly involved in the transition.

In order to overcome the size-consistency problem while ensuring a physically meaningful redistribution of velocity, here we propose a method that adjusts nuclear velocities based on atomic contributions to the electronic transition. These contributions are derived from either the one-electron transition density matrix or the electron density difference between the two involved electronic states.⁵⁵ The one-electron transition density matrix and the electron density difference provide essential insights into the spatial characteristics of electronic transitions.⁵⁶ Since these quantities are readily available in many electronic structure methods, our approach is widely applicable. By applying population analysis to those matrices, we can identify the most relevant atoms involved in the transition and redistribute kinetic energy accordingly. This ensures that velocity rescaling remains size-consistent and physically meaningful while minimizing unnecessary adjustments to unrelated degrees of freedom.

Accordingly, this paper introduces two velocity rescaling methods based on atomic contributions. The first method is weighted per-atom velocity rescaling, where the kinetic energy is rescaled per atom according to its excitation weight derived from the one-electron transition density matrix or the density difference matrix. A correction factor ensures energy conservation while maintaining physically meaningful momentum adjustments. The second method is a simpler approach based on threshold-based velocity rescaling. This scheme adjusts only the velocities of atoms whose contributions exceed a predefined excitation threshold. By applying velocity rescaling solely to those atoms actively participating in the electronic transition, the redistribution of kinetic energy remains localized to the most relevant nuclear degrees of freedom, preserving both physical accuracy and computational efficiency. The new methods are tested on the excited-state dynamics of the fulvene molecule (Figure 1, left), because previous studies have shown that population decay in the S₁ state is strongly influenced by the choice of velocity rescaling method in SH.³⁰,³¹,³³,⁴⁴ Fulvene is therefore a prototype system to validate the methodologies proposed. Additionally, we examine the excited-state dynamics of 1H-1,2,3-triazole (triazole) in order to validate our approach on a system containing heteroatoms (Figure 1, right).

The remainder of the paper is structured as follows. Next, we briefly introduce the population analysis used to identify the atomic contributions to electronic transitions and then explain the two proposed velocity rescaling methods. The excited-state

<figure><img src="image_12.png" /><figcaption>**Figure 1.** Structures of fulvene (left) and 1H-1,2,3-triazole (right). Carbon atoms are shown in gray, hydrogen in white, and nitrogen in blue.</figcaption></figure>

dynamics for fulvene and triazole are presented, focusing on the impact of the velocity rescaling algorithm on electronic and nuclear dynamics. Finally, we summarize our findings.

## 2. METHODOLOGY

### 2.1. Atomic Contributions to Electronic Transitions.

We assess how individual atoms contribute to an electronic transition by assigning atomic charges associated with the redistribution of electron density. To this end, we use two complementary approaches.

#### 2.1.1. Approach 1: Electron and Hole Decomposition.
In the first approach, we decompose the transition into electron and hole components, which can be analyzed independently. This decomposition can be performed on either the one-electron transition density matrix (1TDM) or the density difference (DD) matrix. For clarity, we denote both generically as T. Both matrices are valid choices, and we have evaluated both in our analysis.

Within multiconfigurational frameworks—consistent with the methodology employed in this study—the electron and hole transition density matrices ($D^e$ and $D^h$) are defined as

$$D_{pq}^{e} = \sum_{r} T_{rp} T_{rq}, \quad D_{pq}^{h} = \sum_{r} T_{pr} T_{qr} \qquad (2)$$

where p, q, r index the active space orbitals.

To obtain atomic contributions, we perform a Löwdin population analysis on both $D^e$ and $D^h$. The Löwdin scheme is chosen for its orthogonalization properties, which improve interpretability and reduce (though not eliminate) basis set dependence.⁵⁷ This analysis yields atomic electron and hole charges, $q^e$ and $q^h$.

The total contribution of an atom A to the transition is computed as

$$q_A = |q_A^e| + |q_A^h| \qquad (3)$$

The use of absolute values ensures that the electron and hole contributions do not cancel, thus reflecting the full extent of atomic participation.

We refer to this method as qeh. If based on the 1TDM, it is labeled qeh/1TDM; if based on the DD matrix, it is labeled qeh/DD.

#### 2.1.2. Approach 2: Transition Charge Analysis.
In the second approach, we apply a Löwdin analysis directly to the 1TDM to obtain transition charges, $q^{tr.}$.⁵⁶,⁵⁸ Each atom's contribution to the excitation is then:

$$q_A = |q_A^{tr}| \qquad (4)$$
This method is referred to as *qtr*.

Although we exclusively use Löwdin charges in this work, other definitions are also feasible. For example, one could compute the electrostatic potential from a given density matrix (e.g., the 1TDM) and fit atomic charges to reproduce it. Exploring the effects of these alternative schemes is left for future work.

## 2.2. Excitation-Weighted Velocity Rescaling ($v_w$).
This method redistributes the available kinetic energy after a surface hop according to each atom's contribution to the electronic transition.

We first define a normalized weight for each atom A:

$$c_A = \frac{q_A}{\sqrt{\sum_B q_B^2}} \qquad (5)$$

where $q_A$ is the atomic contribution from either eq 3 or eq 4.

The total kinetic energy available to the system undergoing an electronic transition is

$$E_{\text{kin}}^{\text{avail}} = \frac{1}{2} \sum_{A} w_{A} M_{A} \mathbf{v}_{A}^{2} \qquad (6)$$

where $w_A = c_A^2$, $M_A$ is the atomic mass of atom A, and $v_A$ is the velocity of atom A.

To redistribute this energy, we define a per-atom velocity scaling factor:

$$f_A = \sqrt{\frac{E_{\text{kin}}^{\text{new}}}{E_{\text{kin}}^{\text{old}}}} \times [1 + \alpha(w_A - \text{mean}(w))] \qquad (7)$$

Here, $\alpha$ is a tuning parameter controlling the strength of the weighting. If $\alpha = 0$, the scaling is uniform and corresponds to eq 1. Larger $\alpha$ values emphasize atoms with higher contributions.

However, directly scaling velocities by $f_A$ alone does not conserve total kinetic energy. To ensure energy conservation, we first compute the kinetic energy associated with the scaled velocities:

$$E_{\text{kin}}^{*} = \frac{1}{2} \sum_{A} M_{A} (f_{A} \mathbf{v}_{A})^{2} \qquad (8)$$

We then apply a global normalization factor:

$$N = \sqrt{\frac{E_{\text{kin}}^{\text{new}}}{E_{\text{kin}}^{*}}} \qquad (9)$$

(note that the kinetic energy after the hop is calculated as $E_{\text{kin}}^{\text{new}} = E_{\text{tot}} - E_{\text{pot}}^{\text{new}}$). The final rescaled velocities are given by

$$\mathbf{v}_{A}^{\text{new}} = \mathbf{v}_{A}^{\text{old}} \times f_{A} \times N \qquad (10)$$

## 2.3. Excitation-Thresholded Velocity Rescaling ($v_t$).

This method restricts velocity rescaling to atoms with significant contributions to the excitation. An atom is considered significant if its contribution satisfies:

$$q_A > n \times \max(\mathbf{q}) \qquad (11)$$

where $n \in [0, 1]$ is a user-defined threshold and $q_A$ is the atomic contribution from either eq 3 or eq 4.

After a surface hop, velocity rescaling is applied only to the subset of atoms that meet this criterion, using the standard rescaling procedure (see eq 1). This targeted redistribution

ensures that kinetic energy adjustments are limited to the atoms most involved in the electronic transition.

## 2.4. Summary of the Methods.
In total, we have defined three approaches to determine atomic contributions to a given electronic transition: one based on electron-hole charges derived from the 1TDM ($qeh/1TDM$), another using electron-hole charges from the DD matrix ($qeh/DD$), and a third based on transition charges ($qtr$).

Additionally, we have proposed two approaches for velocity rescaling after a hop: the excitation-weighted rescaling ($v_w$) and the excitation-thresholded rescaling ($v_t$). Each of these rescaling methods can be combined with any of the three charge analysis approaches, resulting in six possible method variants, see Table 1. For example, applying excitation-thresholded rescaling with the electron-hole charges computed from the 1TDM is denoted as $v_w(qeh/1TDM)$.

<table><caption>Table 1. Summary of the Methods Resulting from the Combination of Two Different Velocity Rescaling Schemes—$v_w$ and $v_t$—with Three Approaches to the Determination of Atomic Contributions to Electronic Transitions: $qeh/1TDM$, $qeh/DD$, and $qtr$</caption><thead><tr><th></th><th>qeh/1TDM</th><th>qeh/DD</th><th>qtr</th></tr></thead><tbody><tr><th>v<sub>w</sub></th><td>v<sub>w</sub>(qeh/1TDM)</td><td>v<sub>w</sub>(qeh/DD)</td><td>v<sub>w</sub>(qtr)</td></tr><tr><th>v<sub>t</sub></th><td>v<sub>t</sub>(qeh, 1TDM)</td><td>v<sub>t</sub>(qeh, DD)</td><td>v<sub>t</sub>(qtr)</td></tr></tbody></table>

# 3. VALIDATION SYSTEMS

Our approach is validated on two representative systems: fulvene and 1H-1,2,3-triazole (Figure 1). The choice of fulvene is motivated by the previous work of Ibele and Curchod⁴⁴ who reported significant differences in the nonadiabatic dynamics when using SH, depending on whether the velocity adjustment was made along the velocity direction (eq 1) or the NACV direction. These observations were further supported by Toldo et al.³³ who explored additional algorithms for velocity rescaling. Fulvene's small size and ultrafast excited-state relaxation dynamics make it an ideal benchmark for testing novel approaches to nonadiabatic dynamics.³³,⁴⁴,⁵⁹ The second system is 1H-1,2,3-triazole, a five-membered aromatic heterocycle containing three nitrogen atoms. This molecule is known for its rich photochemistry and ultrafast nonradiative decay pathways.⁶⁰⁻⁶² Previous studies have shown that its $S_1$ excited state is dissociative, leading to elongation and rupture of the N2-N3 bond. This process relaxes the system to the ground state, potentially releasing molecular nitrogen ($N_2$) and forming products such as ethanimine.⁶⁰,⁶¹

We emphasize that the primary objective of our study is not to provide a comprehensive mapping of the photochemical behavior or product distribution of 1H-1,2,3-triazole, but rather to evaluate the performance of our velocity rescaling algorithm within the first 150 fs of excited-state dynamics.

Although the two approaches proposed for velocity rescaling ($v_w$ and $v_t$) can be particularly suited for larger systems, here we prioritized assessing its performance and generalizability in small molecules like fulvene and triazole.

Specifically for fulvene, we tested three variants of the excitation-weighted rescaling method (Section 2.2): $v_w(qeh/1TDM)$, using $\alpha = 0.25$, 0.50, and 1.00 (see eq 7 for the definition of the $\alpha$ parameter); $v_w(qeh/DD)$, using $\alpha = 0.50$ and 1.00; and $v_w(qtr)$, using $\alpha = 0.50$ and 1.00. We also considered one variation of the excitation-thresholded method,
$v_i(qeh, 1TDM)$, which applies a threshold of $n = 0.30$ (see eq 11). In addition, we benchmarked our approaches ($v_w$ and $v_i$) against two commonly used methods: velocity adjustment along the NACV direction (d) and the standard velocity rescaling method described in eq 1 ($v_{full}$). We further compared our simulations with results from the ab initio multiple spawning (AIMS) method,⁶³⁻⁶⁵ which provides a robust reference due to its insensitivity to velocity rescaling schemes.

For 1H-1,2,3-triazole, we applied the same three variants of the excitation-weighted rescaling method — $v_w(qeh/1TDM)$, $v_w(qeh/DD)$, and $v_w(qtr)$ — using $\alpha = 0.50$ and $\alpha = 1.00$ for each approach, and compared the results with both the NACV-based velocity adjustment (d) and the standard rescaling method ($v_{full}$).

## 4. COMPUTATIONAL DETAILS

**4.1. Electronic Structure.** The electronic properties required to propagate the nuclear dynamics of fulvene—including energies, nuclear gradients, and NACVs—were computed at the SA(3)-CASSCF(6,6)/6-31G* level of theory. The active space comprised six electrons distributed over six molecular orbitals, encompassing the full $\pi$ system: three pairs of $\pi$ and $\pi*$ orbitals (see Figure S1).

For triazole, the calculations were performed at the SA(7)-CASSCF(10,8)/6-31G* level, with an active space of ten electrons in eight orbitals: two nonbonding ($n$) orbitals, two pairs of $\pi$ and $\pi*$ orbitals, and one pair of $\sigma$ and $\sigma*$ orbitals localized on the N2-N3 bond (see Figure S2). Initially, the SA(5)-CASSCF(10,8)/6-31G* level was tested for triazole. However, many trajectories failed early in the simulations due to convergence issues. To improve stability, we increased the number of averaged states to SA(7), which resolved these problems and allowed the simulations to proceed reliably.

All electronic structure calculations were performed using the OpenMolcas program, version 23.10 for SH dynamics and version 24.10 for AIMS dynamics.⁶⁶

Population analyses were carried out using the internal WFA module in OpenMolcas. This module directly yields the transition charges (eq 4), as well as the hole and electron charges (eq 3), computed from both the 1TDM and the DD matrix.

**4.2. Nonadiabatic Dynamics.** Initial coordinates and velocities were stochastically sampled from a Wigner distribution corresponding to uncoupled harmonic oscillators, constructed from a harmonic frequency analysis at the ground state optimized geometry. For fulvene, 200 initial conditions—comprising geometries, velocities, and initial electronic states—were selected within an excitation energy window of 3.7–4.3 eV, targeting excitation to the $S_1$ state (see absorption spectrum in Figure S5). The nuclear dynamics were propagated with a time step of 0.2 fs for a total simulation time of 50 fs. For triazole, 119 initial conditions were selected based on an excitation energy range of 5.75–6.25 eV, also corresponding to excitation to the $S_1$ state (see Figure S6). The nuclear dynamics were propagated with a time step of 0.2 fs for a total simulation time of 150 fs.

For each system, all SH simulations were launched from the same set of geometries and velocities. To ensure stochastic consistency across simulations, different velocity rescaling schemes were applied using identical random seeds. The fewest switches algorithm was used in the SH simulations, incorporating the local diabatization approach⁵³,⁵⁴ for prop-

agating the electronic coefficients (even when NACVs were employed for velocity adjustment) with the Granucci-Persico energy-based decoherence correction¹⁷ using the empirical parameter of 0.1 au. In cases of frustrated hops, the momentum direction was preserved. The SH simulations were propagated using the SHARC molecular dynamics package.⁶⁷

The AIMS dynamics on fulvene were performed for a subset of 170 initial conditions (geometries, velocities, and initial electronic states) drawn from the SH ensemble. The quantum amplitudes were propagated with a full diagonal propagator in the adiabatic basis with norm-preserving interpolation time-derivative couplings.⁶⁸ A time step of 10 atomic time units (approximately 0.24 fs) was used and the absolute value of the time-derivative coupling above 0.0039 au⁻¹ was used as a spawning criterion. The velocity rescaling of newly spawned trajectories was performed for the full velocity vector, as it was observed in previous studies⁴⁴ that this has no effect on AIMS dynamics. The width of the Gaussians was chosen to be 4.7 and 22.7 1/bohr² for H and C, respectively. The AIMS simulations were run with the Pyspawn code.⁶⁹⁻⁷¹

## 5. RESULTS AND DISCUSSION

**5.1. Fulvene.** The photodynamics of fulvene have been extensively studied in the literature.³³,⁴⁴,⁴⁸,⁵⁹,⁷²⁻⁷⁶ Upon photoexcitation to the $S_1$ state, it undergoes ultrafast excited-state decay through two primary nonadiabatic pathways mediated by conical intersections with the $S_0$ state. The dominant pathway is driven by the stretching of the C1–C2 bond (see Figure 1), leading to a strongly sloped conical intersection. The nuclear wavepacket passes through this intersection, undergoes reflection, and recrosses, resulting in a stepwise decay of the $S_1$ population. The alternative pathway involves a twist of the same C1–C2 bond, where a fully twisted geometry corresponds to a peaked conical intersection, while a partial twist results in a smoother, less pronounced slope. An extended $S_1/S_0$ seam enables population transfer across a range of geometries, from planar to highly twisted. Simulations and theoretical studies consistently indicate that most of the population decays via the sloped conical intersection, while a smaller fraction follows the twisting mechanism. Therefore, in terms of nuclear dynamics, the most important internal coordinate to investigate is the C1–C2 bond length, followed by the two cis dihedrals (C3–C2–C1–H2 and C6–C2–C2–H1) and the two trans dihedrals (C3–C2–C1–H1 and C6–C2–C2–H2). In Figure 2, we present the $S_1$ state population decay calculated using SH with different velocity adjustments after a hopping event, as well as for AIMS. Specifically, we compare velocity-direction adjustment using the full kinetic energy available ($v_{full}$), NACV-direction adjustment (d), the three different variants of excitation-weighted velocity rescaling ($v_w$), and excitation-thresholded velocity rescaling ($v_i(qeh, 1TDM)$).

Upon excitation to the $S_1$, the energy of the $S_1$ and $S_0$ states approaches each other, reaching a conical intersection after 7 fs, which effectively transfers about 80% of the population from $S_1$ to $S_0$. However, after approximately 12 fs, a back transfer of population to the $S_1$ state is observed. The coupling region is encountered again after 24 fs. At around 24 fs, the population decays again to the $S_0$ state, and after approximately 32 fs, a second repopulation of the $S_1$ state is observed, though less pronounced this time.
<figure><img src="image_13.png" /><figcaption>**Figure 2.** (a) Time-resolved S₁ state population of fulvene over the first 50 fs, computed using surface hopping with different velocity adjustment algorithms: velocity-direction adjustment using the full kinetic energy available (*v*<sub>full</sub>), NACV-direction adjustment (*d*), three different flavors of the excitation-weighted velocity rescaling (*v*<sub>w</sub>), and excitation-thresholded velocity rescaling (*v*<sub>t</sub>(qeh, 1TDM)). The S₁ population, obtained with ab initio multiple spawning (AIMS) is also showed. Shaded regions indicate 95% confidence intervals (Γ) for the NACV-direction adjustment results, calculated as Γ = *p* ± 1.96 × √[(*p*(1−*p*))/*N*<sub>traj</sub>], where *p* is the state population and *N*<sub>traj</sub> = 197 is the number of trajectories. (b): Magnified view of the S₁ population across all approaches.</figcaption></figure>

All the simulations exhibit a similar initial population decay. The recrossing and back transfer to S₁ occur in all approaches, but the extent of population transfer depends strongly on the nonadiabatic dynamics method employed. While AIMS predicts that only about 24% of the population remains in S₁ between 16 and 24 fs, the repopulation process in SH occurs slightly later, leading to higher excited-state population plateaus. The final S₁ population reaches approximately 31% for the **d** adjustment, 29% for the **v**<sub>w</sub>(qeh/1TDM) adjustment, 31% for the **v**<sub>w</sub>(qeh/DD) adjustment, 27% for the **v**<sub>w</sub>(qtr) adjustment, 40% for the **v**<sub>t</sub>(qeh, 1TDM) adjustment, and 39% for the **v**<sub>full</sub> adjustment.

By using the **d** adjustment as a reference for SH dynamics, we observe that both **v**<sub>full</sub> and **v**<sub>t</sub>(qeh, 1TDM) significantly overestimate the repopulation to the S<sub>1</sub> state, as previously noted in refs 33 and 44. The **v**<sub>w</sub> approaches underestimate the repopulation in comparison to **d** adjustment; however, **v**<sub>w</sub>(qeh/1TDM) and **v**<sub>w</sub>(qeh/DD) show much better agreement with the **d** dynamics, representing a notable improvement in the population dynamics. It is gratifying to see that for all **v**<sub>w</sub> approaches, the S<sub>1</sub> population over the entire time interval of

the dynamics falls within the 95% confidence interval (shaded region in Figure 2) of the **d** dynamics, indicating consistency between these methods within statistical uncertainty.

In addition, we have analyzed the internal consistency between the classical populations (Figure S7) and the quantum populations (Figure 2), and from this comparison, one can observe that internal consistency is preserved. Note that the classical populations are computed as Nᵢ(t)/N<sub>traj</sub> where Nᵢ(t) is the number of trajectories in which the active state at time t is i, and N<sub>traj</sub> is the total number of trajectories. In contrast, the quantum populations are calculated as the average over all trajectories of the state probabilities |Cᵢ(t)|², where Cᵢ(t) are the electronic adiabatic coefficients.

We use AIMS dynamics as a benchmark, as it has been previously shown that it is insensitive to rescaling schemes.<sup>44,75</sup> One of the reasons for this observation is the standard spawning algorithm<sup>65,77</sup> as implemented in PySpawn, as it ensures that the new trajectory is created at the point of maximum coupling, which generally is expected to coincide with the point of minimal energy difference which also minimizes the extend of kinetic energy rescaling. The difference between the SH and AIMS population dynamics is particularly evident during the first decay, where the AIMS S<sub>1</sub> population drops to 20% after 12 fs. The subsequent repopulation up to 25% in AIMS is reproduced in the SH simulations that use **d** or **v**<sub>w</sub> approaches. The AIMS population trace at this point does not lie within the 95% confidence interval of the **d** dynamics; however, this discrepancy is primarily due to differences in the initial deactivation behavior (within the first 12 fs). It is important to note that the first deactivation is minimally, if at all, affected by the choice of rescaling scheme after a hop or spawn as the rescaling only affects the behavior of the trajectory after the hop and does not alter the initial hopping probability to energetically lower lying states. Therefore, the observed deviations can be attributed mainly to differences in the number of initial conditions, the degree of convergence achieved, and the fundamental conceptual distinctions between the SH and AIMS methods. Despite these differences, the excellent agreement in the repopulation phase between the **v**<sub>w</sub> dynamics and AIMS serves as strong validation for the velocity rescaling scheme.

**Table 2** presents the total number of hops (S<sub>1</sub> → S<sub>0</sub> transitions), back-hops (S<sub>0</sub> → S<sub>1</sub> transitions), and frustrated hops across all SH simulations. Additionally, we report the total number of trajectories considered in our analysis. Trajectories where total energy conservation was violated (total energy variation was superior to |0.20| eV) or where issues in the electronic structure calculations occurred were excluded. To provide a normalized comparison, values in parentheses represent the number of hops, back-hops, and frustrated hops events per trajectory, e.g., $\frac{\text{Number of hops}}{\text{Number of trajectories}}$.

The results clearly indicate that the increased repopulation of the S<sub>1</sub> state observed with **v**<sub>full</sub> and **v**<sub>t</sub>(qeh, 1TDM) is primarily due to an excessive number of back-hopping events. Specifically, **v**<sub>full</sub> and **v**<sub>t</sub>(qeh, 1TDM) exhibit 0.410 and 0.426 back-hops per trajectory, respectively—substantially higher than those observed for the other methods, which do not exceed 0.229. This suggests that these approaches provide excess kinetic energy for nonadiabatic transitions.

In contrast, the **d** adjustment, as well as the excitation-weighted velocity rescaling (*v*<sub>w</sub>) methods, significantly reduce the number of back-hops. The reduction in back-hopping for
<table><caption>**Table 2.** Number of Hops (N<sub>hop</sub>), Back-Hops (N<sub>back hop</sub>), Frustrated Hops (N<sub>frust</sub>), and Total Trajectories Considered (N<sub>traj</sub>) for Fulvene under Different Velocity Adjustment Schemes Following a Hopping Event<sup>ab</sup></caption><thead><tr><th></th><th>N<sub>hop</sub></th><th>N<sub>back hop</sub></th><th>N<sub>frust</sub></th><th>N<sub>traj</sub></th></tr></thead><tbody><tr><td><b>d</b></td><td>226 (1.147)</td><td>39 (0.198)</td><td>52 (0.264)</td><td>197</td></tr><tr><td><b>v<sub>w</sub><sup>full</sup></b></td><td>260 (1.300)</td><td>82 (0.410)</td><td>0 (0.000)</td><td>200</td></tr><tr><td><b>v<sub>w</sub><sup>qeh/1TDM</sup>, α = 1.00)</b></td><td>225 (1.148)</td><td>41 (0.209)</td><td>53 (0.270)</td><td>196</td></tr><tr><td><b>v<sub>w</sub><sup>qeh/1TDM</sup>, α = 0.50)</b></td><td>226 (1.141)</td><td>39 (0.197)</td><td>53 (0.268)</td><td>198</td></tr><tr><td><b>v<sub>w</sub><sup>qeh/1TDM</sup>, α = 0.25)</b></td><td>228 (1.146)</td><td>42 (0.211)</td><td>61 (0.306)</td><td>199</td></tr><tr><td><b>v<sub>w</sub><sup>qtr</sup>, α = 1.00)</b></td><td>226 (1.147)</td><td>37 (0.188)</td><td>48 (0.244)</td><td>197</td></tr><tr><td><b>v<sub>w</sub><sup>qtr</sup>, α = 1.00)</b></td><td>230 (1.156)</td><td>41 (0.206)</td><td>43 (0.216)</td><td>199</td></tr><tr><td><b>v<sub>w</sub><sup>qeh/DD</sup>, α = 1.00)</b></td><td>235 (1.187)</td><td>47 (0.237)</td><td>48 (0.242)</td><td>198</td></tr><tr><td><b>v<sub>w</sub><sup>qeh/DD</sup>, α = 0.50)</b></td><td>232 (1.184)</td><td>45 (0.229)</td><td>47 (0.240)</td><td>196</td></tr><tr><td><b>v<sub>t</sub><sup>qeh</sup>, 1TDM)</b></td><td>257 (1.318)</td><td>83 (0.426)</td><td>1 (0.005)</td><td>195</td></tr></tbody></table>

<sup>a</sup>The velocity adjustments include velocity-direction adjustment using the full kinetic energy reservoir (v<sub>full</sub>), adjustment along the nonadiabatic coupling vectors (d), various flavors of excitation-weighted velocity rescaling (v<sub>w</sub>), and excitation-thresholded velocity rescaling (v<sub>t</sub>(qeh, 1TDM)). <sup>b</sup>Values in parentheses indicate the number of hopping, back hopping, and frustrated hopping events per trajectory (e.g., N<sub>hop</sub>/N<sub>traj</sub>).

these approaches can be attributed to the lower available kinetic energy for hopping, as constrained by their respective rescaling mechanisms. The derivation of the kinetic energy available when using NACV to adjust velocities is provided in Appendix A (see eq 22). Overall, our results suggest that the d and v<sub>w</sub> approaches provide a more balanced description of nonadiabatic transitions by reducing spurious back-hoppings. A key observation is that the excitation-weighted velocity rescaling (v<sub>w</sub>) induces a hopping behavior that closely resembles the d approach. The v<sub>w</sub> methods not only maintain a lower number of back-hops but also introduce a controlled number of frustrated hops, preventing excessive transitions and thus better agreement with physically realistic nonadiabatic dynamics.

To compare the kinetic energy available for hopping, we extracted the back-hopping geometries and velocities from simulations using the d adjustment. We then computed the available kinetic energy using the full kinetic energy reservoir (v<sub>full</sub>), adjustment along the nonadiabatic coupling vectors (d), various flavors of excitation-weighted velocity rescaling (v<sub>w</sub>), and excitation-thresholded velocity rescaling (v<sub>t</sub>(qeh, 1TDM)). The results are presented in Figure 3, where each x-component corresponds to a different geometry. We analyzed a total of 38 geometries.

As we can see, the kinetic energy available for hopping in the v<sub>w</sub> approach is often underestimated compared to that obtained using the NACV-direction adjustment (d). However, it shows better agreement than the other approaches. For both d and v<sub>w</sub>, the kinetic energy available for hopping remains below 1.0 eV (except for one geometry in the d method). In contrast, when considering the full kinetic energy reservoir (v<sub>full</sub>), the available kinetic energy is typically above 1.5 eV. When focusing only on atoms that contribute the most to the electronic transition to compute the kinetic energy (v<sub>t</sub>(qeh, 1TDM)), we observe a significant improvement; however, in many cases, the values are still underestimated compared to the d adjustment.

<figure><img src="image_14.png" /><figcaption>**Figure 3.** Kinetic energy available for hopping across different approaches for fulvene. Each x-component corresponds to a back-hopping geometry and velocity extracted from simulations using the NACV-direction adjustment. The kinetic energy is computed using four different methods: d (eq 22), v<sub>w</sub> (eq 6), the full kinetic energy reservoir v<sub>full</sub>, and v<sub>t</sub>, which considers only the atoms contributing most to the electronic transition when computing the kinetic energy.</figcaption></figure>

Although the results obtained for the electronic dynamics are quite satisfactory, the first proposed method (v<sub>w</sub>) prompted us to investigate whether it also captures and improves the nuclear dynamics compared to, for instance, v<sub>full</sub>. To this end, in Figure 4 we plot the time evolution of the average bond length between atoms C2 and C1 of fulvene (recall Figure 1), obtained from SH with different velocity adjustment algorithms. Additionally, Figure 4b shows the difference in bond length relative to the d simulations. We observe that the C1–C2 stretching oscillates more slowly when using the v<sub>full</sub> method compared to the others. This behavior is expected, as v<sub>full</sub> redistributes the total energy equally across all atoms, without considering their individual contributions to the nonadiabatic process. In contrast, the NACV-direction adjustment and excitation-weighted velocity rescaling methods (v<sub>w</sub>) account for the varying atomic contributions, leading to a more localized energy redistribution. Since the C1–C2 stretching mode is known to actively participate in the S<sub>1</sub>/S<sub>0</sub> transition, a stronger local adjustment enhances the kinetic energy in this specific coordinate, resulting in faster oscillations of the C1–C2 bond.

Importantly, the v<sub>w</sub> also appears to enhance nuclear dynamics, as it preferentially distributes energy to atoms actively involved in electronic transitions while minimizing unnecessary energy redistribution to those that are not. This seems to improve the size-consistency problem associated with the SH approach using v<sub>full</sub>, without requiring explicit NACV calculations.

Next, we analyzed the normalized contribution of the NACV between the S<sub>1</sub>/S<sub>0</sub> states per atom for various back-hopping geometries extracted from d adjustment simulations, using the same geometries as in Figure 3. We then computed the average of these coefficients across all geometries and compared them to the average coefficients obtained from eq 5 for the different v<sub>w</sub> approaches.

Figure 5 compares the average normalized contribution of the NACV per atom for the different back-hopping geometries extracted from simulations using the d adjustment with the average coefficients calculated according to eq 5 for the various v<sub>w</sub> approaches. The x-axis represents the atoms (see Figure 1).
<figure><img src="image_15.png" /><figcaption>**Figure 4.** (a) Time evolution of the average bond length between atoms C2 and C1 of fulvene (see Figure 1), computed using surface hopping with different velocity adjustment algorithms: velocity-direction adjustment using the full kinetic energy reservoir ($v_{\text{full}}$), NACV-direction adjustment (**d**), and excitation-weighted velocity rescaling ($v_w$). (b) Difference in bond length relative to the **d** simulations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 5.** Comparison of the atomic contribution in the velocity adjustment for the **d** and $v_w$ in fulvene's dynamics. We present the average normalized contribution of the NACV per atom for different back-hopping geometries extracted from simulations using the **d** adjustment is compared with the average coefficients calculated according to eq 5 for the different $v_w$ approaches. The x-axis represents individual atoms (see Figure 1). The y-axis shows the average coefficient values and the normalized nonadiabatic coupling contribution for the 39 back-hopping geometries obtained through **d** calculations.</figcaption></figure>

The y-axis shows the average of the coefficients (eq 5) and the normalized nonadiabatic coupling contribution from the **d**

dynamics for the 38 back-hopping geometries obtained through **d** calculations. Our objective is to assess whether the atomic contributions used to adjust the velocity in $v_w$ align with those in the **d**. From Figure 5, we observe that the individual atomic contribution $v_w(qeh/DD)$ aligns best with the atomic contribution obtained from **d**, and in general, the qualitative agreement is satisfactory.

As mentioned previously, the dominant pathway is driven by the stretching of the C1–C2 bond. For $v_w(qeh/DD)$, the contributions from C1 and C2 are slightly stronger than those from **d**. This may explain why the stretching bond oscillations shown in Figure 4 are somewhat stronger in $v_w(qeh/DD)$ compared to **d**. For $v_w(qeh/1\text{TDM})$, the contribution from C1 is overstated, while the contribution from C2, which has the highest contribution in **d**, is understated. This discrepancy may lead to a cancellation error, but the oscillation of the C1–C2 stretching shows good agreement with **d**.

Using $v_w(qtr)$, the contributions from C1 and C2 are both understated. As shown in Figure 4, the oscillations of this bond exhibit a slight delay, particularly at longer times. This effect is even more pronounced with the $v_{\text{full}}$ approach, where the oscillations are substantially delayed compared to other methods due to the equal distribution of kinetic energy among all atoms.

In summary, regarding electronic dynamics, the $v_w$ method, in particular the $v_w(qeh/DD)$, demonstrates substantial improvement compared to the other approaches, aligning closely with the NACV-direction adjustment. Additionally, we observe that the electronic dynamics remain robust with respect to the chosen $\alpha$ value. However, the electronic dynamics obtained with the $v_t$ method is quite similar to that observed with $v_{\text{full}}$, especially because the system still retains excess kinetic energy, resulting in an excessive number of back-hoppings. We note that the threshold-based velocity rescaling procedure is inherently system-dependent. In small molecules such as fulvene, where the kinetic energy contributions from heavy atoms (e.g., carbon) are relatively uniform and the influence of light atoms (e.g., hydrogen) is non-negligible, the method offers limited improvement. As a result, varying the “n” parameter has little to no impact on the outcome. In contrast, for larger chromophores or multichromophoric systems, this approach can be more effective, as it enables selective weighting of specific fragments or atomic types.

**5.2. H-1,2,3-Triazole.** The temporal evolution of the $S_0$, $S_1$, and $S_2$ state populations, comparing different velocity-direction adjustment methods, is shown in Figure 6a. Specifically, the full kinetic energy adjustment ($v_{\text{full}}$), the NACV-direction adjustment (**d**), and three variants of excitation-weighted velocity rescaling ($v_w(qeh/1\text{TDM})$) for $\alpha = 1.00$ are shown. Results for $\alpha = 0.50$ can be found in Figure S3.

The $S_1$ and $S_0$ population curves exhibit an inverse trend: an initial plateau is observed during the first 30 fs of the simulation, which can be attributed to skeleton relaxation along the N2–N3 bond length alternation. Subsequently, the $S_1$ population decays to the ground state, while a small fraction is transferred to the $S_2$ state in all approaches. The magnified view of the $S_2$ population (Figure 6b) reveals that the $v_{\text{full}}$ method results in a slightly higher $S_2$ population. This suggests that providing more kinetic energy facilitates back-hopping to the $S_2$ state, compared to the other approaches.

The $S_0$ and $S_1$ populations over the entire time interval of the dynamics show that all approaches fall within the 95% confidence interval (shaded region) of the **d** dynamics,
<figure><img src="image_17.png" /><figcaption>**Figure 6.** (a) Time-resolved state populations of 1H-1,2,3-triazole over the first 150 fs, computed using different velocity adjustment algorithms: velocity-direction adjustment with full kinetic energy available ($v_{full}$), NACV-direction adjustment (d), and three variations of excitation-weighted velocity rescaling ($v_w$). The $S_0$ state is represented by solid lines (—), $S_1$ by dashed lines (-), and $S_2$ by dotted lines (···). Shaded regions indicate 95% confidence intervals ($\Gamma$) for the NACV-direction adjustment results, calculated as $\Gamma = p \pm 1.96 \times \sqrt{\frac{p(1-p)}{N_{\text{traj}}}}$, where p is the state population and $N_{\text{traj}} = 110$ is the number of trajectories. (b) Magnified view of the $S_2$ population across all approaches.</figcaption></figure>

indicating consistency between these methods within statistical uncertainty. However, the $S_2$ population for the $v_{full}$ approach, for a significant portion of the dynamics, does not fall within the 95% confidence interval compared to the d simulations.

Analogous to fulvene, Table 3 presents the total number of hopping events ($S_1 \rightarrow S_0$, $S_2 \rightarrow S_0$, and $S_2 \rightarrow S_1$ transitions), back-hopping events ($S_0 \rightarrow S_1$, $S_0 \rightarrow S_2$, and $S_1 \rightarrow S_2$ transitions), and frustrated hopping events across all SH simulations in H-1,2,3-triazole. Additionally, we report the total number of trajectories analyzed, excluding those where total energy conservation was violated (total energy variation was superior to |0.50| eV) or where electronic structure calculations encountered issues. For a normalized comparison, values in parentheses represent the number of hopping, back-hopping, and frustrated hopping events per trajectory.

The results indicate that the increased population of the $S_2$ state observed with $v_{full}$ primarily arises from a higher number of back-hopping events from $S_1 \rightarrow S_2$ compared to the other approaches. Notably, $v_{full}$ is the only approach where an $S_0 \rightarrow S_2$ back-hopping event was observed. All $v_w$ simulations showed improvements in the number of hopping and back-hopping events, yielding results closer to those obtained with the NACV-based approach (d) compared to $v_{full}$. In particular, the $v_w(qeh/1TDM)$ method produced results that most closely aligned with those of d.

To evaluate the kinetic energy available for hopping, we extracted $S_0 \rightarrow S_1$ back-hopping geometries and velocities from simulations using the d adjustment. As shown in Figure 7, the

<figure><img src="image_18.png" /><figcaption>**Figure 7.** Kinetic energy available for hopping across different approaches for triazole. Each x-component corresponds to an $S_0 \rightarrow S_1$ back-hopping geometry extracted from simulations using the NACV-direction adjustment. The kinetic energy is computed using four methods: d, $v_w$, $v_{traj}$, and $v_{full}$.</figcaption></figure>

<table><caption>**Table 3.** Total Number of Hops ($N_{\text{hop}}$), Back-Hops ($N_{\text{back hop}}$), and Frustrated Hops ($N_{\text{frust}}$) Events, along with the Total Number of Trajectories Considered ($N_{\text{traj}}$), for Triazole under Different Velocity Adjustment Schemes<sup>a</sup></caption><thead><tr><th></th><th>N<sup>S1</sup><sub>hop</sub>→S<sub>0</sub></th><th>N<sup>S0</sup><sub>backhop</sub></th><th>N<sup>S2</sup><sub>hop</sub>→S<sub>1</sub></th><th>N<sup>S0</sup><sub>backhop</sub></th><th>N<sup>S2</sup><sub>hop</sub>→S<sub>1</sub></th><th>N<sup>S1</sup><sub>backhop</sub>→S<sub>2</sub></th><th>N<sub>frust</sub></th><th>N<sub>traj</sub></th></tr></thead><tbody><tr><td>d</td><td>75 (0.682)</td><td>21 (0.191)</td><td>0 (0.000)</td><td>0 (0.000)</td><td>3 (0.027)</td><td>4 (0.036)</td><td>29 (0.264)</td><td>110</td></tr><tr><td>v<sub>full</sub></td><td>81 (0.750)</td><td>33 (0.305)</td><td>1 (0.009)</td><td>1 (0.009)</td><td>15 (0.139)</td><td>17 (0.157)</td><td>1 (0.009)</td><td>108</td></tr><tr><td>v<sub>w</sub>(qeh/1TDM, α = 1.00)</td><td>73 (0.676)</td><td>22 (0.204)</td><td>0 (0.000)</td><td>0 (0.000)</td><td>11 (0.102)</td><td>12 (0.111)</td><td>15 (0.139)</td><td>108</td></tr><tr><td>v<sub>w</sub>(qeh/1TDM, α = 0.50)</td><td>72 (0.679)</td><td>21 (0.198)</td><td>1 (0.009)</td><td>0 (0.000)</td><td>8 (0.075)</td><td>9 (0.085)</td><td>12 (0.113)</td><td>106</td></tr><tr><td>v<sub>w</sub>(qtr, α = 1.00)</td><td>80 (0.748)</td><td>28 (0.262)</td><td>0 (0.000)</td><td>0 (0.000)</td><td>6 (0.056)</td><td>7 (0.065)</td><td>8 (0.075)</td><td>107</td></tr><tr><td>v<sub>w</sub>(qtr, α = 1.00)</td><td>72 (0.692)</td><td>25 (0.240)</td><td>0 (0.000)</td><td>0 (0.000)</td><td>8 (0.077)</td><td>8 (0.077)</td><td>11 (0.106)</td><td>104</td></tr><tr><td>v<sub>w</sub>(qeh/DD, α = 1.00)</td><td>72 (0.686)</td><td>24 (0.229)</td><td>1 (0.009)</td><td>0 (0.000)</td><td>9 (0.086)</td><td>10 (0.095)</td><td>11 (0.105)</td><td>105</td></tr><tr><td>v<sub>w</sub>(qeh/DD, α = 0.50)</td><td>72 (0.692)</td><td>23 (0.221)</td><td>0 (0.000)</td><td>0 (0.000)</td><td>8 (0.077)</td><td>9 (0.086)</td><td>11 (0.106)</td><td>104</td></tr></tbody></table>

<footnote><sup>a</sup>Values in parentheses indicate the number of events per trajectory (e.g., $\frac{N_{\text{hop}}^{S1}-N_{\text{hop}}^{S0}}{N_{\text{traj}}}$).</footnote>
$v_w$ approach generally overestimates the available kinetic energy compared to NACV-based calculations, though occasional underestimations occur. Nevertheless, $v_w$ demonstrates better overall agreement than the other methods. While we did not perform dynamics using the $v_t(qeh, 1TDM)$ approach, we estimated its available kinetic energy for comparison. For both $d$ and $v_w$ the available kinetic energy remains below 1.0 eV, except for two geometries in the $d$ method. In contrast, the full kinetic energy reservoir ($v_{full}$) consistently yields values above 2 eV.

Finally, Figure 8 compares the average normalized atomic contribution of the NACV per atom for various $S_0 \to S_1$ back-

<figure><img src="image_19.png" /><figcaption>**Figure 8.** Comparison of the atomic contributions in velocity adjustment for NACV ($d$) and $v_w$ in triazole dynamics. We present the average normalized contribution of the NACV per atom for different back-hopping geometries extracted from simulations using the $d$ adjustment, compared with the average coefficients calculated according to eq 5 for the different $v_w$ approaches. The $x$-axis represents individual atoms, where atoms (see Figure 1). The $y$-axis shows the average coefficient values and the normalized NACV contribution for the 31 $S_0 \to S_1$ back-hopping geometries obtained through $d$ calculations.</figcaption></figure>

hopping geometries extracted from $d$ adjustment simulations. The $x$-axis represents individual atoms. The $y$-axis shows the average coefficient values from eq 5 and the normalized NACV contribution for the 31 analyzed geometries. As seen in Figure 8, the $v_w(qeh/DD)$ approach most closely reproduces the atomic contributions obtained from $d$. Notably, this was also observed in fulvene.

In summary, the investigated velocity adjustment methods show similar population dynamics over time for triazole, with minimal impact on the nuclear dynamics as well (see Section S3). However, the $v_w$ method displayed a number of hops, back-hops, and frustrated hops that align more closely with the $d$ compared to $v_{full}$. Additionally, the kinetic energy available for hopping in $v_w$ also shows much better agreement with $d$ compared to $v_{full}$.

## 6. CONCLUSIONS

In this study, we addressed the size-consistency issue in the surface hopping method by proposing two novel velocity rescaling approaches. Both methods adjust nuclear velocities based on atomic contributions to electronic transitions, which are derived from either the one-electron transition density matrix or the density difference matrix. The first method, *excitation-weighted velocity rescaling*, redistributes kinetic energy

across atoms based on their individual contributions, using a normalized coefficient derived from the population analysis of the one-electron transition density matrix or density difference matrix. The second method, *excitation-thresholded velocity rescaling*, adjusts velocities only for atoms that exceed a predefined threshold.

The methods were tested to simulate the excited-state dynamics of fulvene and 1H-1,2,3-triazole. The results were compared against those resulting from the standard velocity-direction adjustment, where kinetic energy is redistributed equally across all atoms and all atoms contribute to the kinetic energy available for hopping, as well as from the nonadiabatic coupling vector-direction adjustment. Our results show that *excitation-weighted velocity rescaling* produces results comparable to nonadiabatic coupling vector-direction adjustment. This approach reduces back-hopping events and provides a more balanced description of nonadiabatic transitions. The population dynamics for both fulvene and 1H-1,2,3-triazole show that all variants of the *excitation-weighted velocity rescaling* method fall within the 95% confidence interval of the nonadiabatic coupling vector-direction adjustment dynamics, indicating consistency between these methods within statistical uncertainty. Additionally, it improves the description of nuclear dynamics by allocating more kinetic energy to atoms actively involved in nonadiabatic transitions.

The *excitation-thresholded velocity rescaling* method was tested on fulvene and yielded results similar to the velocity-direction adjustment. This is primarily because the system retains excess kinetic energy, leading to an excessive number of back-hopping events. Although this approach may not be superior for small organic molecules like fulvene, it may be worth exploring for larger chromophores or multichromophoric systems.

Overall, our study demonstrates that considering atomic contributions to electronic transitions provides a viable alternative for velocity adjustment in surface hopping, especially when nonadiabatic coupling vectors are unavailable or too computationally expensive to calculate.

## ■ APPENDIX

### Velocity adjustment along the nonadiabatic coupling direction

The total classical energy must be conserved during hopping between states. If the system hops from an old state to a new state, we must have

$$E_{\text{kin}}^{\text{old}} + E_{\text{pot}}^{\text{old}} = E_{\text{kin}}^{\text{new}} + E_{\text{pot}}^{\text{new}} \quad (12)$$

where $E_{\text{kin}}$ represents the kinetic energy and $E_{\text{pot}}$ represents the potential energy. The potential energy gap between the two states is defined as

$$\Delta E_{\text{pot}} = E_{\text{pot}}^{\text{new}} - E_{\text{pot}}^{\text{old}} \quad (13)$$

Therefore, from the energy conservation equation, the kinetic energy variation is

$$\Delta E_{\text{kin}} = E_{\text{kin}}^{\text{new}} - E_{\text{kin}}^{\text{old}} = -\Delta E_{\text{pot}} \quad (14)$$

To adjust for the change in kinetic energy, the velocity of nucleus A is rescaled as

$$\mathbf{v}_{A}^{\text{new}} = \mathbf{v}_{A}^{\text{old}} + \gamma \frac{\mathbf{d}_{A}^{\text{old, new}}}{M_{A}} \quad (15)$$
where $\mathbf{d}_A^{\text{old,new}}$ is the velocity rescaling direction, in this case, in the direction of the nonadiabatic coupling vector between the states involved in the transition, and $M_A$ is the nuclear mass. Substituting this into the kinetic energy equation (eq 14) leads to

$$a\gamma^2 + b\gamma + \Delta E_{\text{pot}} = 0 \qquad (16)$$

where

$$a = \frac{1}{2} \sum_A \frac{1}{M_A} (\mathbf{d}_A^{\text{old,new}} \cdot \mathbf{d}_A^{\text{old,new}}) \qquad (17)$$

$$b = \sum_A \frac{1}{M_A} (\mathbf{v}_A^{\text{old}} \cdot \mathbf{d}_A^{\text{old,new}}) \qquad (18)$$

being the discriminant of this equation:

$$\Delta = b^2 - 4a\Delta E_{\text{pot}} \qquad (19)$$

Therefore, if $\Delta < 0$, no real solution exists, and the hopping is frustrated. Otherwise, if $\Delta \ge 0$, the possible solutions are

$$\gamma = \begin{cases} \frac{-b + \sqrt{\Delta}}{2a}, & \text{if } -b + \sqrt{\Delta} < -b - \sqrt{\Delta} \\ \frac{-b - \sqrt{\Delta}}{2a}, & \text{otherwise} \end{cases} \qquad (20)$$

Therefore, hopping is only allowed if:

$$\frac{1}{2 \sum_A \frac{(\mathbf{d}_A^{\text{old,new}})^2}{M_A}} \left( \sum_A \mathbf{v}_A^{\text{old}} \cdot \mathbf{d}_A^{\text{old,new}} \right)^2 \ge \Delta E_{\text{pot}} \qquad (21)$$

From the conditions above, the kinetic energy available for hopping is

$$\Delta E_{\text{kin}} \le \frac{1}{2 \sum_A \frac{(\mathbf{d}_A^{\text{old,new}})^2}{M_A}} \left( \sum_A \mathbf{v}_A^{\text{old}} \cdot \mathbf{d}_A^{\text{old,new}} \right)^2 \qquad (22)$$

## ■ ASSOCIATED CONTENT

### ■ Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.5c00737.

Details on the active space used for the studied systems, additional results on the evolution of electronic populations over time, time evolution of some geometric parameters over time for 1H-1,2,3-triazole, and the simulated absorption spectra (PDF)

## ■ AUTHOR INFORMATION

### Corresponding Authors

Eduarda Sangiogo-Gil – University of Vienna, Institute of Theoretical Chemistry, Vienna A-1090, Austria;

● orcid.org/0000-0002-0945-2906;

Email: eduarda.sangiogo.gil@univie.ac.at

Leticia González – University of Vienna, Institute of Theoretical Chemistry, Vienna A-1090, Austria; Research Platform on Accelerating Photoreaction Discovery (ViRAPID), University of Vienna, Vienna 1090, Austria;

● orcid.org/0000-0001-5112-794X;

Email: leticia.gonzalez@univie.ac.at

## Authors

Lea M. Ibele – Aix Marseille University, CNRS, ICR, Marseille 13397, France; ● orcid.org/0000-0002-1434-3843

Richard Bleyer – University of Vienna, Institute of Theoretical Chemistry, Vienna A-1090, Austria

Complete contact information is available at: https://pubs.acs.org/10.1021/acs.jctc.5c00737

## Notes

The authors declare no competing financial interest.

## ■ ACKNOWLEDGMENTS

The authors gratefully acknowledge the University of Vienna for its continuous support and the Vienna Scientific Cluster (VSC) for the generous allocation of computational resources. L.M.I. acknowledges the Centre de Calcul Intensif d'Aix-Marseille for granting access to its high-performance computing facilities.

## ■ REFERENCES

(1) Chapman, S. *Advances in Chemical Physics*. John Wiley & Sons, Ltd., 1992; pp. 423–483.

(2) Barbatti, M. Nonadiabatic dynamics with trajectory surface hopping method. *Wiley Interdiscip. Rev.: Comput. Mol. Sci* **2011**, *1*, 620–633.

(3) Persico, M.; Granucci, G. An overview of nonadiabatic dynamics simulations methods, with focus on the direct approach versus the fitting of potential energy surfaces. *Theor. Chem. Acc* **2014**, *133* (9), 1526.

(4) Wang, L.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011–2015. *J. Phys. Chem. Lett* **2016**, *7*, 2100–2112.

(5) Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N. Understanding the Surface Hopping View of Electronic Transitions and Decoherence. *Annu. Rev. Phys. Chem* **2016**, *67*, 387–417.

(6) Mai, S.; Marquetand, P.; González, L. Nonadiabatic dynamics: The SHARC approach. *Wiley Interdiscip. Rev.: Comput. Mol. Sci* **2018**, *8*, No. e1370.

(7) Wang, L.; Qiu, J.; Bai, X.; Xu, J. Surface hopping methods for nonadiabatic dynamics in extended systems. *Wiley Interdiscip. Rev.: Comput. Mol. Sci* **2020**, *10*, No. e1435.

(8) Cigrang, L. L. E.; et al. Roadmap for Molecular Benchmarks in Nonadiabatic Dynamics. *J. Phys. Chem. A* **2025**, *129*, 7023.

(9) Tully, J. C. Molecular dynamics with electronic transitions. *J. Chem. Phys* **1990**, *93*, 1061–1071.

(10) Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum-Classical Dynamics. *Chem. Rev* **2018**, *118*, 7026–7068.

(11) Jain, A.; Sindhu, A. Pedagogical Overview of the Fewest Switches Surface Hopping Method. *ACS Omega* **2022**, *7*, 45810–45824.

(12) Toldo, J. M.; Do Casal, M. T.; Ventura, E.; Do Monte, S. A.; Barbatti, M. Surface hopping modeling of charge and energy transfer in active environments. *Phys. Chem. Chem. Phys* **2023**, *25*, 8293–8316.

(13) Plasser, F.; Gómez, S.; Menger, M. F. S. J.; Mai, S.; González, L. Highly efficient surface hopping dynamics using a linear vibronic coupling model. *Phys. Chem. Chem. Phys* **2019**, *21*, 57–69.

(14) Zobel, J. P.; Heindl, M.; Plasser, F.; Mai, S.; González, L. Surface Hopping Dynamics on Vibronic Coupling Models. *Acc. Chem. Res* **2021**, *54*, 3760–3771.

(15) Farkhutdinova, D.; Polonius, S.; Karrer, P.; Mai, S.; González, L. Parametrization of Linear Vibronic Coupling Models for Degenerate Electronic States. *J. Phys. Chem. A* **2025**, *129*, 2655–2666.
(16) Granucci, G.; Persico, M.; Zoccante, A. Including quantum decoherence in surface hopping. *J. Chem. Phys* **2010**, *133* (13), 134111.

(17) Granucci, G.; Persico, M. Critical appraisal of the fewest switches algorithm for surface hopping. *J. Chem. Phys* **2007**, *126* (13), 134114.

(18) Jasper, A. W.; Truhlar, D. G. Non-Born-Oppenheimer molecular dynamics of Na···FH photodissociation. *J. Chem. Phys* **2007**, *127*, 194306.

(19) Persico, M.; Granucci, G.; Accomasso, D. *Comprehensive Computational Chemistry* 1st ed.: Yáñez, M.; Boyd, R. J. Eds.; Elsevier: Oxford; 2024 pp. 273–292.

(20) Pieroni, C.; Sangiogo Gil, E.; Ibele, L. M.; Persico, M.; Granucci, G.; Agostini, F. Investigating the Photodynamics of trans-Azobenzene with Coupled Trajectories. *J. Chem. Theory Comput* **2024**, *20*, 580–596.

(21) Belyaev, A. K.; Lasser, C.; Trigila, G. Landau–Zener type surface hopping algorithms. *J. Chem. Phys* **2014**, *140*, 224108.

(22) Mannouch, J. R.; Richardson, J. O. A mapping approach to surface hopping. *J. Chem. Phys* **2023**, *158* (10), 104111.

(23) Karplus, M.; Porter, R. N.; Sharma, R. D. Exchange Reactions with Activation Energy. I. Simple Barrier Potential for (H, H2). *J. Chem. Phys* **1965**, *43*, 3259–3287.

(24) Gray, J. C.; Truhlar, D. G.; Clemens, L.; Duff, J. W.; Chapman, F. M.; Morrell, G. O.; Hayes, E. F. Quasiclassical trajectory calculations compared to quantum mechanical reaction probabilities, rate constants, and activation energies for two different potential surfaces for the collinear reaction H2+I → H+HI, including dependence on initial vibrational state. *J. Chem. Phys* **1978**, *69*, 240–252.

(25) Mukherjee, S.; Barbatti, M. A Hessian-Free Method to Prevent Zero-Point Energy Leakage in Classical Trajectories. *J. Chem. Theory Comput* **2022**, *18*, 4109–4116.

(26) Pieroni, C.; Becuzzi, F.; Creatini, L.; Granucci, G.; Persico, M. 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* **2023**, *19*, 2430–2445.

(27) Janoš, J.; Slaviček, P.; Curchod, B. F. E. Selecting Initial Conditions for Trajectory-Based Nonadiabatic Simulations. *Acc. Chem. Res* **2025**, *58*, 261–270.

(28) Tully, J. C. Nonadiabatic molecular dynamics. *Int. J. Quantum Chem* **1991**, *40*, 299–309.

(29) Herman, M. F. Nonadiabatic semiclassical scattering. I. Analysis of generalized surface hopping procedures. *J. Chem. Phys* **1984**, *81*, 754–763.

(30) Jasper, A. W.; Hack, M. D.; Truhlar, D. G. The treatment of classically forbidden electronic transitions in semiclassical trajectory surface hopping calculations. *J. Chem. Phys* **2001**, *115*, 1804–1816.

(31) Jasper, A. W.; Truhlar, D. G. Improved treatment of momentum at classically forbidden electronic transitions in trajectory surface hopping calculations. *Chem. Phys. Lett* **2003**, *369*, 60–67.

(32) 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**, *17*, 3010–3018.

(33) Toldo, J. M.; Mattos, R. S.; Pinheiro, M. J.; Mukherjee, S.; Barbatti, M. Recommendations for Velocity Adjustment in Surface Hopping. *J. Chem. Theory Comput* **2024**, *20*, 614–624.

(34) Sangiogo Gil, E.; Granucci, G.; Persico, M. Surface Hopping Dynamics with the Frenkel Exciton Model in a Semiempirical Framework. *J. Chem. Theory Comput* **2021**, *17*, 7373–7383.

(35) Sangiogo Gil, E.; Persico, M.; Granucci, G. Frenkel exciton photodynamics of self-assembled monolayers of azobiphenyls. *J. Chem. Phys* **2022**, *157* (16), 161101.

(36) Sangiogo Gil, E.; Giustini, A.; Accomasso, D.; Granucci, G. Excitonic Approach for Nonadiabatic Dynamics: Extending Beyond the Frenkel Exciton Model. *J. Chem. Theory Comput* **2024**, *20*, 8437–8449.

(37) Röhr, M. I. S.; Marciniak, H.; Hoche, J.; Schreck, M. H.; Ceymann, H.; Mitric, R.; Lambert, C. Exciton Dynamics from Strong to Weak Coupling Limit Illustrated on a Series of Squaraine Dimers. *J. Phys. Chem. C* **2018**, *122*, 8082–8093.

(38) Ibele, L. M.; Sánchez-Murcia, P. A.; Mai, S.; Nogueira, J. J.; González, L. Excimer Intermediates en Route to Long-Lived Charge-Transfer States in Single-Stranded Adenine DNA as Revealed by Nonadiabatic Dynamics. *J. Phys. Chem. Lett* **2020**, *11*, 7483–7488.

(39) Accomasso, D.; Granucci, G.; Wibowo, M.; Persico, M. Delocalization effects in singlet fission: Comparing models with two and three interacting molecules. *J. Chem. Phys* **2020**, *152* (24), 244125.

(40) Pedraza-González, L.; Accomasso, D.; Cupellini, L.; Granucci, G.; Mennucci, B. Ultrafast excited-state dynamics of Luteins in the major light-harvesting complex LHCII. *Photochem. Photobiol. Sci* **2024**, *23*, 303–314.

(41) Avagliano, D.; Bonfanti, M.; Garavelli, M.; González, L. QM/MM nonadiabatic dynamics: the SHARC/COBRAMM approach. *J. Chem. Theory Comput* **2021**, *17*, 4639–4647.

(42) Truhlar, D. G.; Muckerman, J. T. Bernstein, R. B. *Atom-Molecule Collision Theory: a Guide for the Experimentalist*; Plenum Press: New York, 1979; pp. 505–566.

(43) Pechukas, P. Time-Dependent Semiclassical Scattering Theory. II. Atomic Collisions. *Phys. Rev* **1969**, *181*, 174–185.

(44) Ibele, L. M.; Curchod, B. F. E. A molecular perspective on Tully models for nonadiabatic dynamics. *Phys. Chem. Chem. Phys* **2020**, *22*, 15183–15196.

(45) Coker, D. F.; Xiao, L. Methods for molecular dynamics with nonadiabatic transitions. *J. Chem. Phys* **1995**, *102*, 496–510.

(46) Hack, M. D.; Jasper, A. W.; Volobuev, Y. L.; Schwenke, D. W.; Truhlar, D. G. Quantum Mechanical and Quasiclassical Trajectory Surface Hopping Studies of the Electronically Nonadiabatic Predisociation of the Å State of NaH2. *J. Phys. Chem. A* **1999**, *103*, 6309–6326.

(47) Kapral, R.; Ciccotti, G. Mixed quantum-classical dynamics. *J. Chem. Phys* **1999**, *110*, 8919–8929.

(48) Do Casal, M. T.; Toldo, J. M.; Pinheiro, M., Jr.; Barbatti, M. Fewest switches surface hopping with Baeck-An couplings. *Open Res. Europe* **2021**, *1*, 49.

(49) Araujo, L.; Lasser, C.; Schmidt, B. FSSH-2: Fewest Switches Surface Hopping with Robust Switching Probability. *J. Chem. Theory Comput* **2024**, *20*, 3413–3419.

(50) Shu, Y.; Zhang, L.; Chen, X.; Sun, S.; Huang, Y.; Truhlar, D. G. Nonadiabatic Dynamics Algorithms with Only Potential Energies and Gradients: Curvature-Driven Coherent Switching with Decay of Mixing and Curvature-Driven Trajectory Surface Hopping. *J. Chem. Theory Comput* **2022**, *18*, 1320–1328.

(51) Hammes-Schiffer, S.; Tully, J. C. Proton transfer in solution: Molecular dynamics with quantum transitions. *J. Chem. Phys* **1994**, *101*, 4657–4667.

(52) Pittner, J.; Lischka, H.; Barbatti, M. Optimization of mixed quantum-classical dynamics: Time-derivative coupling terms and selected couplings. *Chem. Phys* **2009**, *356*, 147–152.

(53) Granucci, G.; Persico, M.; Toniolo, A. Direct semiclassical simulation of photochemical processes with semiempirical wave functions. *J. Chem. Phys* **2001**, *114*, 10608–10615.

(54) Plasser, F.; Granucci, G.; Pittner, J.; Barbatti, M.; Persico, M.; Lischka, H. 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* **2012**, *137* (22), 22A514.

(55) Plasser, F.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. I. Formalism. *J. Chem. Phys* **2014**, *141* (2), 024106.

(56) Plasser, F. TheoDORE: A toolbox for a detailed and automated analysis of electronic excited state computations. *J. Chem. Phys* **2020**, *152* (8), 084108.
(57) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: introduction to Advanced Electronic Structure Theory, 1st ed.; Dover Publications: New York, 1989.

(58) Voityuk, A. A. Fragment transition density method to calculate electronic coupling for excitation energy transfer. *J. Chem. Phys* **2014**, *140* (24), 244117.

(59) Avagliano, D.; Lorini, E.; González, L. Sampling effects in quantum mechanical/molecular mechanics trajectory surface hopping non-adiabatic dynamics. *Phil. Trans. R. Soc. A* **2022**, *380*, 20200381.

(60) Burgess, E. M.; Carithers, R.; McCullagh, L. Photochemical decomposition of 1H-1, 2, 3-triazole derivatives. *J. Am. Chem. Soc* **1968**, *90*, 1923-1924.

(61) Sangiogo Gil, E.; de Araújo, B. B.; Gonçalves, P. F. B. CASPT2, CASSCF and non-adiabatic molecular dynamics (NAMD) studies on the low-lying electronic states of 1H-1,2,3-triazole photolysis. *Phys. Chem. Chem. Phys* **2019**, *21*, 25809-25819.

(62) Ghosh, J.; Gajapathy, H.; Jayachandran, A.; Bernstein, E. R.; Bhattacharya, A. On the electronically nonadiabatic decomposition dynamics of furazan and triazole energetic molecules. *J. Chem. Phys* **2019**, *150* (16), 164304.

(63) Ben-Nun, M.; Martínez, T. J. Nonadiabatic molecular dynamics: Validation of the multiple spawning method for a multidimensional problem. *J. Chem. Phys* **1998**, *108*, 7244-7257.

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

(65) Curchod, B. F. E. *Quantum Chemistry and Dynamics of Excited States*. John Wiley & Sons, Ltd., 2020; pp. 435-467.

(66) Li Manni, G.; et al. The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry. *J. Chem. Theory Comput* **2023**, *19*, 6933-6991.

(67) Mai, S., et al. SHARC4.0: *Surface Hopping Including Arbitrary Couplings — Program Package for Non-Adiabatic Dynamics*. 2025 https://sharc-md.org/.

(68) Meek, G. A.; Levine, B. G. Evaluation of the time-derivative coupling for accurate electronic state transition probabilities from numerical simulations. *J. Phys. Chem. Lett* **2014**, *5*, 2351-2356.

(69) Fedorov, D. A.; Seritan, S.; Fales, B. S.; Martínez, T. J.; Levine, B. G. PySpawn: Software for nonadiabatic quantum molecular dynamics. *J. Chem. Theory Comput* **2020**, *16*, 5485-5498.

(70) Ibele, L. M.; Memhood, A.; Levine, B. G.; Avagliano, D. Ab Initio Multiple Spawning nonadiabatic dynamics with different CASPT2 flavors: A fully open-source PySpawn/OpenMolcas interface. *J. Chem. Theory Comput* **2024**, *20*, 8140-8151.

(71) Avagliano, D. Solvent Effects on Nonadiabatic Dynamics: Ab Initio Multiple Spawning Propagated on CASPT2/xTB Potentials. *J. Chem. Theory Comput* **2025**, *21*, 1905-1915.

(72) Mendive-Tapia, D.; Lasorne, B.; Worth, G. A.; Bearpark, M. J.; Robb, M. A. Controlling the mechanism of fulvene S 1/S 0 decay: switching off the stepwise population transfer. *Phys. Chem. Chem. Phys* **2010**, *12*, 15725-15733.

(73) Mendive-Tapia, D.; Lasorne, B.; Worth, G. A.; Robb, M. A.; Bearpark, M. J. Towards converging non-adiabatic direct dynamics calculations using frozen-width variational Gaussian product basis functions. *J. Chem. Phys* **2012**, *137* (22), 22A548.

(74) Bearpark, M. J.; Bernardi, F.; Olivucci, M.; Robb, M. A.; Smith, B. R. Can fulvene S1 decay be controlled? A CASSCF study with MMVB dynamics. *J. Am. Chem. Soc* **1996**, *118*, 5254-5260.

(75) Ibele, L. M.; Sangiogo Gil, E.; Schürger, P.; Agostini, F. A coupled-trajectory approach for decoherence, frustrated hops and internal consistency in surface hopping. *arXiv*, **2024**.

(76) Gómez, S.; Spinlove, E.; Worth, G. Benchmarking non-adiabatic quantum dynamics using the molecular Tully models. *Phys. Chem. Chem. Phys* **2024**, *26*, 1829-1844.

(77) Ibele, L. M.; Nicolson, A.; Curchod, B. F. Excited-state dynamics of molecules with classically driven trajectories and Gaussians. *Mol. Phys* **2020**, *118*, No. e1665199.

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

CAS BIOFINDER DISCOVERY PLATFORM™

BRIDGE
BIOLOGY AND
CHEMISTRY FOR
FASTER
ANSWERS

Analyze target relationships,
compound effects, and disease
pathways

Explore the platform

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