# Ab Initio Multiple Spawning nonadiabatic dynamics with different CASPT2 flavors: A fully open-source PySpanw/OpenMolcas interface

Lea Ibele, Arshad Memhood, Benjamin Levine, Davide Avagliano

► To cite this version:

Lea Ibele, Arshad Memhood, Benjamin Levine, Davide Avagliano. Ab Initio Multiple Spawning nonadiabatic dynamics with different CASPT2 flavors: A fully open-source PySpanw/OpenMolcas interface. Journal of Chemical Theory and Computation, 2024, 20 (18), pp.8140-8151. 10.1021/acs.jctc.4c00855 . hal-05288272

HAL Id: hal-05288272

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

Submitted on 29 Sep 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>
Ab Initio Multiple Spawning nonadiabatic
dynamics with different CASPT2 flavors: A
fully open-source PySpann/OpenMolcas
interface

Lea M. Ibele,\*,† Arshad Memhood,‡ Benjamin G. Levine,‡ and Davide Avagliano\*,¶

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

‡Institute for Advanced Computational Science and Department of Chemistry, 11794, Stony
Brook, NY, USA

¶Chimie ParisTech, PSL University, CNRS, Institute of Chemistry for Life and Health Sciences (iCLEHS UMR 8060), 75005 Paris, France

E-mail: lea-maria.ibele@universite-paris-saclay.fr; davide.avagliano@chimieparistech.psl.eu

Abstract

We introduce an interface between PySpann, a simulation package to run ab initio
multiple spawning (AIMS) nonadiabatic dynamics, and OpenMolcas, a software package
to perform multiconfigurational perturbations theory (CASPT2) electronic structure
calculations. Our interface allows us to exploit all the functionalities of the two codes:
the modular and efficient Python implementation of the AIMS algorithm and the
extensive analysis tools offered by PySpann, with the cutting-edge implementation
of CASPT2 equations in OpenMolcas, including the recently introduced analytical
gradients and different flavors. Both are fully open-source and free of charge, making
the following implementation unique in the current plethora of software for nonadiabatic
dynamics. This represents an important step towards a wider application of AIMS-based
nonadiabatic dynamics combined with high-accuracy excited-state calculations. The
importance and the need for such an implementation are demonstrated by application
to the ultrafast relaxation of fulvene from S<sub>1</sub> to S<sub>0</sub>, which is drastically affected by the
potential energy surface on which the nuclear wavepacket is propagated. Additionally,
the decay is influenced by the CASPT2 flavor adopted, posing interesting questions in
the choice of one over the other and opening the door to deeper studies on the effect of
CASPT2 formulations in nonadiabatic dynamics.

Introduction

The advancement of nonadiabatic dynamics simulations relies on the parallel development of their two essential components: the way the nuclei are propagated and the electronic structure method used.¹,² These developments must also face the well-known and unavoidable accuracy vs. cost trade-off that leads to choosing one method over another.³,⁴ Regarding the nuclear propagation among different potential energy surfaces, the space of methods is spanned by trajectory-based mixed quantum-classical methods, such as trajectory surface hopping, on the one side³, and fully quantum wavefunction propagation, such as multiconfigurational time-dependent Hartree, on the other extreme of the spectrum. A family of methods that can be placed in the middle are Gaussian-based approaches. They rely on the expansion of the nuclear wavefunction in a basis of Gaussian functions leading to a trajectory-like approach while retaining a nuclear wavefunction in the formalism. The propagation of the trajectory basis functions (TBFs), whether they follow quantum, mixed quantum-classical, or classical trajectories, differentiates the three most popular Gaussian-based methods, variational multiconfigurational Gaussian,⁵,⁶ multiconfigurational Ehrenfest⁷,⁸ and full and Ab initio multiple spawning,⁹,¹⁰ respectively. Despite their different formalisms, in the limit of a sufficiently large Gaussian basis and an exact evaluation of the matrix elements, all
these methods will converge to an exact solution. Ab initio multiple spawning (AIMS)¹¹,¹²
has risen to a popular choice of nonadiabatic dynamics methods for the description of
the photochemistry of medium-sized molecules. It has successfully been used to model
the excited state dynamics of important biomolecules,¹³⁻¹⁵ and to directly simulate time-
resolved spectra.¹⁶⁻²¹ In particular, thanks to its extensions to the description of intersystem
crossing,²²,²³ external fields,²⁴⁻²⁶ tunnelling²⁷ and complex environments,²⁸,²⁹ AIMS has
become a method of wide applicability that offers a good balance between accuracy and cost.
Depending on the system and implementation of the AIMS algorithm, it can provide accurate
dynamics at a cost comparable to or only slightly higher than that of independent trajectory
methods.³⁰ The largest advantage of AIMS over independent trajectory methods is it being
based on controlled approximations and not needing any sort of ad-hoc corrections to correctly
capture i.e. decoherence effects.³,³¹ The original most widely used implementation of AIMS,
FMS90, has been integrated within the Molpro program package³² as well as TeraChem.
Additionally, interfaces have been created for Columbus, GAMESS, and Mopac. There has
been much recent interest in the development of broadly available software packages for
nonadiabatic molecules dynamics simulations,³³⁻³⁹ following the example of the electronic
structure community. In this spirit, the Levine group has developed PySpawn,³⁶ an open-
source software package for AIMS, which is designed to be highly extensible for integration
with various electronic structure codes. Developed in Python, this program offers an intuitive
interface for both software developers and end-users. It includes an interactive analysis module
designed to facilitate the broader adoption of AIMS. Moreover, PySpawn employs a novel
task-based implementation of the AIMS algorithm and utilizes modern structured output
formats, enabling efficient deployment on shared high-performance computing resources.

From the electronic structure point of view, second-order multiconfigurational perturbation
theory applied to reference wavefunction (CASPT2),⁴⁰,⁴¹ generated at the Complete Active
Space Space Self-Consistent Field (CASSCF),⁴² has been shown to push the limit of the
accuracy with a computational cost still affordable to obtain accurate dynamics in the
ultrafast time scale for relatively small chromophores.43–48 The computational bottleneck
is given by the choice of the active space (AS) and its size, which influences the accuracy
of the calculation; alas often the Sword of Damocles is on the computational feasibility of
nonadiabatic simulations using such a method to calculate electronic energies on-the-fly
at each time step. However, the continuous developments of the community have been
partially overcoming some of the limitations to the applicability of CASPT2 in dynamics. In
particular, the recently introduced analytical gradients49–53 and extended multi-state variation
of CASPT2,54–56 specifically developed to deal with strongly correlated wavefunctions, have
now made simulations possible that were previously fantastical. One of the main software to
run and advance multiconfigurational perturbation theory is, without doubt, OpenMolcas,57,58
due to its extensive and continued implementation made by its community, its accuracy,
reliability and its being fully open-source. However, although AIMS and CASPT2 lie among
the highest standards in terms of accuracy/cost ratio, no software including the AIMS
algorithm is currently interfaced with OpenMolcas to exploit all its potential, and all the
CASPT2 flavors implemented, to propagate AIMS nonadiabatic dynamics. This work fills
this gap and introduces a new interface between PySpawn and OpenMolcas to run AIMS
nonadiabatic dynamics at the state-of-the-art CASSCF//CASPT2 level. One of the main
strengths of this interface is offering an open-source implementation of AIMS with a high level
of accuracy in the calculations of electronic energies. The interface allows the exploitation of
all the already existing implementations inside OpenMolcas, including the recently introduced
CASPT2 analytical gradients, which represents a notable speed-up in the excited states
dynamics simulation at such a level of theory. However, the interface was designed flexibly so
it will be easy to include all the future features that will be introduced in OpenMolcas by its
broad, and constantly active community. In the following sections, we will first revise the
PySpawn algorithm and structure, to offer a full overview of how the software propagates
the nuclear wavefunction. Then, we will introduce the main characteristic of the interface
and we will finally show an exemplary application. Thanks to our interface we could test
the effect of different and recent CASPT2 implementations in the nonadiabatic dynamics of fulvene following the excitation to $S_1$, comparing with the well-known CASSCF one. These capabilities make the following implementation unique in the current plethora of software for nonadiabatic dynamics. This represents an important step towards a wider application of AIMS-based nonadiabatic dynamics combined with high-accuracy excited states calculations.

AIMS nonadiabatic dynamics with Py spawn

AIMS is derived from the formally exact full multiple spawning (FMS) framework. Starting from the Born-Huang expansion of the molecular wavefunction

$$ \Psi(\mathbf{r}, \mathbf{R}, t) = \sum_{J} \phi_{J}(\mathbf{r}; \mathbf{R}) \chi^{J}(\mathbf{R}, t) \qquad (1) $$

which writes the wavefunction as a linear combination of electronic eigenstates $\phi_J(\mathbf{r}; \mathbf{R})$ and their time-dependence is grouped into the nuclear expansion coefficients $\chi_J(\mathbf{R}, t)$, FMS proposes to express the nuclear wavefunction in a basis of Gaussian functions.

$$ \chi_J(\mathbf{R}, t) = \sum_k C_k^{(J)}(t) \tilde{\chi}_k^J(\mathbf{R}; \bar{\mathbf{R}}_k^J(t), \bar{\mathbf{P}}_k^J(t), \bar{\alpha}) \qquad (2) $$

where the trajectory basis functions (TBFs) $\tilde{\chi}_k^J(\mathbf{R}; \bar{\mathbf{R}}_k^J(t), \bar{\mathbf{P}}_k^J(t), \bar{\alpha})$ are defined 3N dimensional frozen Gaussians that evolve classically, i.e. their central position and momentum are propagated following Newton's equations of motion. Thus, a de-facto moving grid is created on the support of which the time-dependent Schrödinger equation is solved. Inserting this expansion into the TDSE allows to obtain evolution equations for the time-dependent coefficients as

$$ \dot{\mathbf{C}}^J = -i \mathbf{S}_{JJ}^{-1} \left[ (\mathbf{H}_{JJ} - i\dot{\mathbf{S}}_{JJ}) \mathbf{C}^J + \sum_{I \neq J} \mathbf{H}_{JI} \mathbf{C}^I \right] \qquad (3) $$
where $\mathbf{S}_{JJ}$ is the overlap matrix between TBFs on state $J$ and $\mathbf{H}_{JJ}$ and $\mathbf{H}_{JI}$ are intra- and interstate Hamiltonian matrices that largely govern the population transfer between TBFs on the same or different states, respectively.

Each TBF is associated with an electronic state in which it evolves adiabatically. Nona-diabatic population transfer occurs via the coupling between TBFs on different electronic states mediated by the nonadiabatic couplings. The most characteristic feature of AIMS is the time-dependence of the size of the Gaussian basis. If a TBF encounters a region of high nonadiabaticity, it will spawn a new function onto the corresponding state. This function initially has an amplitude of zero but is fully coupled to the other TBFs in the simulation and ensures a smooth transfer of amplitude between the two coupled electronic states.

The application of FMS to study molecules in their full dimensionality is hindered by the cost of the evaluation of the Hamiltonian matrix elements as they require multidimensional integrals that contain the potential energy surfaces and nonadiabatic coupling terms over the full molecular configuration space. Therefore, the framework of AIMS makes the simulation of molecules accessible by introducing two main approximations. 1) The saddle-point approximation (SPA) alleviates the nuclear-coordinate dependence of the potential energy surfaces or nonadiabatic couplings in the integrand by expanding them in a truncated Taylor series centered at the centroid position of the two TBFs. Usually, a SPA of order zero, SPA0, is used where the truncation is applied after the first term. 2) The independent first generation approximation treats the initial (parent) TBFs describing the nuclear wavefunction (or wavepacket) at time $t = 0$ as uncoupled. Accurate propagation depends on careful treatment of the nonadiabatic couplings, which become singular at conical intersections (CoIns). The first derivative couplings are computed using the norm-preserving interpolation approximation, which ensures accurate propagation even at trivially avoided crossings.⁵⁹ Second derivative couplings are neglected completely. This appears on its face to be an approximation but has actually been shown to accurately account for errors arising from unphysical discontinuities in the adiabatic electronic wave functions.⁶⁰
<figure><img src="image_2.png" /><figcaption>Figure 1: Task-based AIMS algorithm in Py spawn. Each Py Spawn iteration involves updating centroids and marking them for electronic structure calculation if needed, generating a prioritized task queue followed by checking if the maximum simulation time and wall time are reached, executing the highest-priority task, marking TBFs for spawning based on derivative coupling thresholds and creating new TBFs as needed, propagating quantum amplitudes, and generating output for restart and data analysis.</figcaption></figure>
The implementations of AIMS in Py spawn uses a task-based approach (Figure 1), lever-
aging an asynchronous propagation of the classical and quantum variables. This approach
recognizes that AIMS simulations involve three hierarchical processes: creating a time-
dependent basis of TBFs through classical propagation and spawning, performing the required
electronic structure calculations at centroid geometries to build the Hamiltonian, and prop-
agating quantum amplitudes. These processes are organized into the so-called simulation
data structure which enables trajectories to propagate independently. The simulation object
encapsulates all simulation state information, including trajectory objects for each TBF,
centroid objects for inter-TBF calculations, quantum propagation data, and a prioritized task
queue. Each trajectory object contains classical positions, momenta, and electronic structure
results. This structure allows TBFs to propagate independently and asynchronously, enabling
different trajectories to exist at different simulation times and efficient calculation of centroids
and quantum amplitudes. Centroid objects are structured similarly to trajectory objects,
and their computations provide the off-diagonal Hamiltonian elements. These computations
only require corresponding trajectory data, while quantum amplitude propagation depends
on all trajectory and centroid data up to the current time. The task queue is a central object
in Py Spawn, enabling the management of tasks based on their prerequisites and the potential
to advance the quantum amplitudes. This queue can contain tasks corresponding to various
simulation times, provided the prerequisite information is available. Py Spawn defines three
primary tasks: propagating a TBF by one time step, backpropagating a TBF by one time
step, and computing a centroid at a particular time step. Each of these tasks entails a single
electronic structure calculation. The propagate and backpropagate tasks involve updating
the classical parameters of TBFs either forward or backward in time, respectively, using
the velocity Verlet algorithm for integration. These tasks require only the current state of
the individual TBF, independent of other TBFs or the quantum amplitudes. The centroid
computation task requires the positions and momenta of two TBFs at the same time point.
The simulation proceeds by computing the highest-priority task, and spawning new TBFs if
derivative coupling thresholds are exceeded. This ensures that the simulation adapts to new
quantum pathways dynamically. Finally, output for restart and data analysis is generated,
ensuring robust data management. Py spawn 's structure allows for fine-grained restarts and
straightforward parallelization. While the current implementation is not yet parallelized
across tasks, the structure is designed to facilitate this in future releases. Each task, typically
an electronic structure calculation, can be distributed across processors or nodes, which can
be re-assigned a new task upon completion, even if others are still computing their tasks,
enhancing computational efficiency and enabling large-scale simulations. By decoupling the
propagation of classical and quantum variables, tasks may be prioritized based on dependency
and the potential to extend quantum propagation. For the full details on the implementations
and the code, we refer to the original publication of Py Spawn.³⁶

Py Spawn / OpenMolcas interface

Our work benefited from the modular structure of Py spawn, which was designed with the
idea of interfacing more software for electronic structure calculations. The interface with
OpenMolcas is structurally identical to the original interface with TeraChem. Analogously,
the required parameters for the electronic structure calculations must be parsed by the
user when submitting the job, namely method and basis set, active space (AS), number
of roots, multiplicity, and scratch directory. Before submitting the simulation, an initial
guess for the orbitals (an INPORB file) must be provided by the user, to ensure that a
good guess with the proper active space is selected. At each time step, Py spawn parses the
current geometry and the wavefunction at the previous time, associated with the specific
trajectory object. The previous CASSCF wavefunction is used as the initial guess for
the orbitals and either the previous CASSCF or CASPT2 wavefunction will be used to
compute the overlap at two consecutive steps, required for the future calculation of the
time-derivative coupling (TDC). Then Py spawn writes the input file, sets up the environment
for OpenMolcas, and launches the software. OpenMolcas is called each time the electronic structure calculation needs to be performed, i.e. propagation and backpropagation of the TBFs and the centroids. The communications between the two codes are handled by a new internal module in PySpawn, which creates a simulation object containing the tasks to be performed and calls OpenMolcas. All the internal programs of OpenMolcas required are called in a single execution: GATEWAY, SEWARD, RASSCF, CASPT2, ALASKA, and RASSI. Analytical gradients are available at CASSCF and CASPT2 levels thanks to the recent implementation included in OpenMolcas. The interface returns energies, gradients, and wavefunction, and computes the overlap matrices at two consecutive steps. This is calculated with RASSI, and the CASPT2 wavefunction is stored by PySpawn and used again later at the following time step of the trajectory object. The TDC will be later computed using these CASPT2 wavefunction overlaps. Similarly, PySpawn also stores the CASSCF wavefunction to be used again in the following step as the initial guess for the orbitals in the following calculations, either propagation or backpropagation. Saving both CASSCF (JobIph file) and CASPT2 (JobMix file) wavefunctions at each time steps, makes possible post-production analysis on the wavefunctions, i.e. calculating the degree of correction at CASPT2 level over the CASSCF wavefunction or the presence of intruder states, by simply running RASSI on specific geometries extracted by the analysis tool, like representative geometries or when a trajectory crashed.

Both CASSCF and CASPT2 levels are available and can be requested, as well as the different PT2 flavors, by the user at the time of the submission. Only state-average (SA-) CASSCF is possible by default and the available CASPT2 flavors are multistate (MS- CASPT2$^{61}$), extended multi-state (XMS-CASPT2$^{54}$), rotated multi-state (RMS-CASPT2$^{56}$) and extended dynamically weighted (XDW-CASPT2$^{55,62}$). Since its development, MS- CASPT2 represents one of the highest standards of accuracy for a manifold of excited states calculated at the Franck-Condon region of the simulation of static absorption spectra.⁶³ However, its application in nonadiabatic dynamics has been limited by the non-smoothness
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic illustration of the interface between Py spawn and OpenMolcas. The electronic structure software is called when it is needed to propagate and backpropagate the trajectory basis functions and the centroid. OpenMolcas computes energies, gradients, and time-derivative couplings at the centroid between two trajectories. Each electronic structure call runs consequently in a single job all OpenMolcas programs needed to compute energies, analytical gradients, and wavefunction overlap at two different geometries.</figcaption></figure>
of the potential energy surfaces when two states are close in energy and strongly mixed,⁵⁰,⁶⁴ like in the proximity of CoIns and avoided crossing between two surfaces. This has been fixed by the introduction of XMS-CASPT2 which solves this problem, arising due to the overestimation of off-diagonal elements of the Hamiltonian that leads the mixing coefficients of the wavefunctions of the different states to be non-smooth, by defining a Hamiltonian that considers the entire Fock operator.⁶⁴ Very recently, XDW-CASPT2 has been introduced, aiming to combine the accuracy of the MS-CASPT2 energies and the smoothness of XMS-CASPT2 potentials. XDW-CASPT2 interpolates between the two methods according to how separated or close the states are, introducing a dynamic weight in the construction of the Fock operators, following either a state-specific or state-average approach. However, the original implementation XDW-CASPT2 would not be suitable for nonadiabatic dynamics when the electronic energies are calculated on-the-fly, since imposing a symmetry is required to avoid the mixing of states of different symmetry. However, in dynamics, no symmetry is imposed. In a more recent formulation,⁶² new ways to define such exponents for the weights are been defined, which also include the symmetry of the states and we thought it would be worth it to test this in nonadiabatic dynamics. In our interface, we left the way and the factors necessary to compute the weights as default and recommended in OpenMolcas to avoid symmetry-related problems, but the user can easily modify the setup. Further slight modifications led to the development of RMS-CASPT2, which improves the description of XDW-CASPT2 as it does not mix states that belong to different irreducible representations, and offers a parameter-free derivation, removing the need to define an empirical parameter to control how to weight between state-average and specific. For all the flavors, in our implementation imaginary and IPEA shifts can be defined in the submission of the dynamics, although values of 0.2 and 0.0 respectively are assigned by default.⁶⁵ As our interface is designed to exploit all the functionality of OpenMolcas, restricted active space, RASSCF, and RASPT2 calculations are also theoretically possible, to include a bigger active space.⁶⁶ However, due to the large number of single-point calculations needed to propagate long enough
dynamics and compute all the quantities required by PySpawn, an accurate cost/benefit
analysis is recommended before running such simulations, so this is the reason why by default
CASSCF/CASPT2 calculations is setup. The interface is included in PySpwawn, which is
available on Git at https://github.com/blevine37/pySpawn17.

Fulvene S<sub>1</sub>/S<sub>0</sub> relaxation with different CASPT2 flavors

The photophysics of fulvene has been extensively studied due to the interesting shape of the
potential energy surfaces (PESs) of the ground and the first excited singlet state.⁶⁷⁻⁷⁰ The
two surfaces cross twice, with one sloped CoIn associated with the stretch of the C=CH₂
bond and one peaked CoIn following the motion of the dihedral angle, still involving the
same moiety. The interesting features are connected to the different relaxation profiles to the
ground state, highly characterized at the CASSCF level, with different nonadiabatic methods.
When the wavepacket meets the sloped CoIn and decays to the ground state, it can hit the
CoIn again, and part of the population gets reflected to S₁. For this feature, fulvene was
recently proposed as a molecular Tully model for nonadiabatic dynamics.³⁰ Alternatively,
encountering the peaked CoIn, the wavepacket relaxes to the ground state without reflection.
In the last years, several methods and parameters have been tested on fulvene,⁷¹⁻⁷⁴ which
represents the perfect playground also to test the PESs at different CASPT2 flavors. Indeed,
no dynamics have been run along different CASPT2 surfaces and extensively compared. We
first ran AIMS/SA-2-CASSCF(6,6)/6-31G* dynamics on a group of 18 initial conditions (IC)
offered to the community to test their method on.³⁰ All the electronic structure calculations
used an active space of 6 electrons and 6 orbitals, corresponding to the π system of the
molecule and the corresponding π* orbitals, and include two singlet electronic states. The
full list of parameters for the dynamics is reported in the Supporting Information. These ICs
differ for the initial geometry and they all have a zero momentum associated. The dynamics
are shown in Figure3, and it is compared with the same level of theory, but with the difference
of the electronic structure software (OpenMolcas vs. MOLPRO) and the implementation of
the AIMS algorithm (Py spawn vs FMS90). The CASSCF relaxation profile of S₁ is analogous
in the two simulations, although a slightly different amount of the reflected population ( 11
vs 16%), which can be attributed to the different implementations of the AIMS algorithm,
in particular, the way the TDC is computed,⁵⁹ the spawning threshold is defined, the use
of an adaptive time step and the momentum rescaled after a spawn. The difference in
population between the two nonadiabatic dynamics codes and electronic structure codes is
not surprising as very small differences in gradients or couplings can have a large impact
on the nonadiabatic dynamics. We also note that the number of initial conditions, 18, is
somewhat limited despite the faster convergence of AIMS simulations in comparison with
conventional independent trajectory approaches.⁷¹,⁷⁵ As each initial condition creates a branch
of multiple interacting TBFs, over longer times the simulations are expected to converge well
as at the end of the propagation the number has grown to a total of 125 TBFs. However, we
expect smaller quantitative differences between the two CASSCF dynamics, in particular
at early times, when the number of TBFs is still smaller. We are certain that increasing
the number of initial conditions would lead to a full convergence between the results, but
given the good qualitative agreement, we are not concerned with the ICs set. Additionally,
we compared the presently computed population profile with a profile computed via the
original Py Spawn/TeraChem interface (Figure S2 of the Supporting Information), showing
an equivalent decay obtained with the two electronic structure software. When propagating
on the support of energies, gradients, and couplings calculated on-the-fly with any flavor of
CASPT2, the ultrafast decays within the first 40fs, using the same ICs, are clearly distinct
with respect to the CASSCF one. Three main differences can be noticed: i) the initial decay
occurs slightly later in time, after 10 fs with respect to the 8 of CASSCF, ii) the initial
relaxation transfers only a small portion of the population to the ground state and with
different extents according to the formulation used (maximum 35% of XDW-CASPT2), iii)
only RMS-CASPT2 shows some degree of reflection of population, while the other dynamics
<figure><img src="image_4.png" /><figcaption>Figure 3: Population of the first excited state of fulvene in the first 40fs of dynamics. SA-2-CASSCF(6,6)/6-31g* profile with our interface and the one reported in ref$^{30}$ with FMS90/MOLPRO implementation are in purple are compared with the profile obtained at different CASPT2 flavors.</figcaption></figure>

show a constant population of S₁ after the initial decay. For completeness, the population profile along the first 100fs at the CASPT2 levels is shown in the Supporting Information (Figure S3). We limited to reporting XMS-CASPT2 and RMS-CASPT2 profiles, as they are the only two dynamics that ensure excellent energy conservation for all the TBFs along the 100fs, while, in some of the TBFs propagated along MS-CASPT2 and XDW-CASPT2 surfaces, jumps in total energy can be observed after passing the CoIns. It can be observed that the main deactivation of the population with all CASPT2 flavors occurs step-wise over longer times. The stronger population transfer to the ground state happens after around 75 fs. At this point, both XMS-CASPT2 and RMS-CASPT2 show a reflection of around 10% of the population comparable to the population trace of CASSCF. After 100 fs, both CASPT2 flavors predict more than 70% of S₁ population having decayed to the ground state. We note that smaller number of 18 initial conditions does not guarantee a fully converged dynamics, hindering a purely quantitative comparison between the different dynamics. However, it is still sufficient to see the stark changes when switching to CASPT2 dynamics and particularly
allows an interesting comparison, as the identical initial conditions and AIMS implementation
are used, highlighting the differences arising purely from the change of electronic structure
method. We wanted to investigate more in-depth the differences between the CASPT2 and
CASSCF PESs and their coupling, leading to such different time evolution. For this reason,
we first computed the static scans along the two coordinates leading to the two degeneracy
points (Figure 4). We linearly interpolated between the ground state optimized geometry
and the CoIns optimized at CASSCF level.⁷⁰ The two coordinates sampled in the two scans
are the ones that represent the coordinate promoting the reflection (CoIn<sub>sloped</sub>, Figure 4a)
and the full decay (CoIn<sub>peaked</sub>, Figure 4b) to the ground state. The sloped CoIn is mainly
characterized by an elongation of the C=CH<sub>2</sub> bond, while the peaked CoIn occurs at a 90°
torsion of the CH<sub>2</sub> group. However, it should be stressed that both these CoIns are minimum
energy CoIns and might not be representative of their full seams of CoIns. Calculating the
overlap between the CASSCF and CASPT2 wavefunctions along the scans confirms that the
order of states matches for each of the points. CASSCF overestimates S<sub>1</sub> energies at the
Franck-Condon point along both coordinates and the wavepacket, initialized with the same
momentum, would thus approach much faster the intersections found at a similar energies
as with CASPT2. Already by looking at the scans, we could anticipate different dynamics
along the CASPT2 surfaces. Along the reaction coordinate leading to the sloped CI, we can
observe all the CASPT2 minima of S<sub>1</sub> far from the intersection between the two surfaces,
with respect to the minimum at the CASSCF level. This could lead to a slower deactivation,
as well as the very sloped profile after the CoIns, could indicate a strong reflection after the
wavepacket crosses it. In the other reaction coordinate, we cannot observe a crossing point
located around the CASSCF CoIn, at any of the CASPT2 levels, which could lead to closing
this deactivation channel at this level of theory.

To understand in more detail the dynamics around the two CoIns, we picked from the
ensemble of initial conditions used to run AIMS/SA-2-CASSCF dynamics one that hits
the sloped CoIns and gives the highest reflection from S₀ to S₁ (IC15) and one the decays
<figure><img src="image_5.png" /><figcaption>Figure 4: Scan along the coordinates connecting the Franck-Condon region and the two conical intersections leading to the relaxation to the ground state.</figcaption></figure>
completely to the ground state after evolving around the peaked CoIn (IC8). In both cases, the dynamics are initiated with zero momentum. CASSCF dynamics show quite different ultrafast dynamics around the two CoIns with respect to the CASPT2 ones. When starting from IC8 (Figure 5C), the wavepacket evolving on the CASSCF surface spawns after 8 fs the first time, and a few fs after a second time, transferring completely population to the ground state. When the wavepacket propagates on any of the CASPT2 profiles, it does not encounter the peaked CoIn, never entering a spawning region in the first 40 fs. By looking at the relative energies between S₀ and S₁ (Figure 6a), we can confirm how the two surfaces cross twice along the CASSCF profile, while they never intersect in any of the CASPT2 flavors. Consequently, the magnitude of the TDC is proximal to zero in these dynamics, while quite high in the CASSCF, promoting the generation of the child trajectories (Figure 6b). Propagating IC15 (Figure 5D), different CASPT2 flavors led to different dynamics, although an overall trend and common difference from the CASSCF one can be seen. In the latter, the wavepacket evolving on S₁ spawns twice in a few fs, analogously to IC8, but the first child spawns at 14 fs inducing a transfer of population back to S₁ of almost 40% of the total population. However, in this IC, the various flavors of CASPT2 actually play a role. The main features of the overall dynamics are in common: slower first spawn with respect to CASSCF, a certain amount of population transferred to the ground state, part of it later reflected to the first singlet state, and finally no further decay to S₀ in the first 40fs. However, slightly different profiles can be noticed. First, although less drastic, the spawn of the first child is slightly scattered from 10.9 to 12.0 fs. More interestingly, both the amount of population decayed to the ground state and the reflected one, vary in the four CASPT2 dynamics, due to the distinct surfaces and couplings between TBFs. XMS-CASPT2 and XDW-CASPT2 show a very limited amount of population decayed to the ground state (around 30%) followed by a small degree of reflection. Much more population is initially transferred at MS-CASPT2 level (ca. 70%), while a full relaxation is observed in the RMS-CASPT2 dynamics. These last two show a comparable S₁ relaxation after reflection. Another important difference worth
<figure><img src="image_6.png" /><figcaption>Figure 5: S₁ population profile for two representative dynamics initialized at two initial conditions that evolve around two distinct conical intersections. Around the peaked CoIn (C), the CASSCF dynamic decays completely to the ground state, while all the CASPT2 dynamics stay on the surface of the first singlet state. After the wavepackets hit the sloped CoIn (D), the population is reflected to S₁, with the extent that differs according to the level of theory employed to calculate electronic energies. All the dynamics are run with (6,6) active space and 6-31g* basis set. The top panel represents the total energy of the parent trajectory basis function along the dynamics around the peaked (A) and sloped (B) conical intersections.</figcaption></figure>
it to be underlined is the different way the population is transferred back to $S_1$. While for CASSCF the child spawns back when it encounters again the CoIn after a few fs, this happens even faster in the even more sloped CoIn in CASPT2 the overlap between the two TBFs is high enough to prevent a second child from being spawned, but not to the population to be transferred back to $S_1$ due to the strong coupling between the two Gaussians. The overall behavior can be better understood by looking again at relative energies and TDC between the two singlet states (Figure6c,d). At CASPT2 levels the two surfaces cross only once, when the CASSCF ones had already crossed twice. It is important to stress that RMS-CASPT2

<figure><img src="image_7.png" /><figcaption>Figure 6: Potential energy profiles for the dynamics evolving around the sloped conical intersection (A) and the peaked one (C), and time derivative coupling calculated along the dynamics (row b: sloped; row d: peaked). Full lines represent $S_0$, $S_1$ energies, and TDC of the parents, while dashed lines are the respective quantities of the spawned child.</figcaption></figure>

and XMS-CASPT2 are the only two flavors that show a perfectly conserved total energy of
the parent TBF in IC15 (Figure5B), while MS-CASPT2 and XDW-CASPT2, with a not drastic but still interesting jump, show some degree of discontinuity once are they passing through the CoIn. On the contrary, along IC8 the total energy is perfectly smooth at any CASPT2 level, as the deactivation channel through the peaked CoIn is closed and indeed the parent TBF does not encounter it and does not manifest any discontinuity in the total energy (Figure5A).

Additionally, we want to assess the influence of the different electronic structure methods on the geometrical evolution of the molecule. In particular, it is interesting to see in which regions of configuration space initially the coupling between electronic states becomes large, driving the deactivation of the molecule. Therefore, Figure 7 shows the C=CH₂ bond length elongation and HCCC torsion angle at which the first spawn occurs in each IC. We note that we plot all first spawns, in the CASPT2 dynamics including those that occur after 43 fs. It has to be noted that the coordinates indicated for the optimized CoIns are obtained at the CASSCF level and as shown in Fig. 4 these do not correspond to CoIns at any CASPT2 level. As discussed in Ref,³⁰ the dynamics are initialized with zero momentum for each initial condition, which is expected to initially lead to rather a stretching of the C=CH₂ bond than the torsion of the CH₂ group. This is reflected in the spawning geometries that overall do not show very large dihedral angles. A general trend can be observed that more spawns occur with all CASPT2 flavors at higher C=CH₂ stretches and also at larger torsion angles. In particular, all CASSCF spawns happen and bond lengths between 1.5 and 1.6Å, while in all CASPT2 dynamics spawns occur outside this range, at smaller as well as longer bond lengths. Curiously, despite only a small fraction of trajectories showing a reflective behavior a majority seems to spawn in the vicinity of the sloped CoIn. This leads to the assumption that while the minimum point of the seam of CoIn dominated by the C=CH₂ bond elongation shows a strongly sloped topology, the same cannot be said for other CoIns along this seam.

We also look closer at the time evolution of these two coordinates along the two initial conditions chosen as representative for the two deactivation pathways above. In figure 8, we
<figure><img src="image_8.png" /><figcaption>Figure 7: C=CH₂ bond length and HCCC torsion angle of the geometries at which the first spawn occurs in all the dynamics. The respective coordinates of the sloped and peaked minimum energy CoIn are indicated by the triangles. We note that the dihedral angles are calculated between 0 and 360°</figcaption></figure>

show the evolution of the C=CH₂ bond length and torsion angle over time for the entire branch of TBFs created from the two initial conditions with CASSCF and all flavors of CASPT2. Each newly spawned TBF is added at the point when it starts being forward propagated. Looking first at the initial condition accessing the peaked CoIn, in the CASSCF dynamics, an initial elongation of the C=CH₂ bond past 1.5Å leads to an initial spawning event closely followed by a second spawn. A majority of the population is almost instantaneously transferred to the first child TBF, at later times a back-and-forth transfer between the two TBFs evolving on the S₀ can be observed. All TBFs show an oscillation of the C=CH₂ bond between almost 1.2 and 1.6Å. Surprisingly, the initial dynamics do not show a strong torsional motion, the HCCC dihedral angle stays mostly constant up until the second spawn, and significant population transfer has occurred.

All the dynamics with different CASPT2 formulations show a similar, if weaker initial elongation of the C=CH₂ bond, oscillating between 1.4 and 1.6Å, showing a similar path in configuration space as explored by the original TBF in CASSCF. As no spawns occur,
only the evolution of the initial TBF can be assessed for this TBF. Curiously, while the evolution of the bond length seems qualitatively the same between all CASPT2 dynamics, the XDW-CASPT2 initial condition remains the slight torsion angle throughout the dynamics, while all other flavors exhibit the fully planar configuration after 43 fs.

The CASSCF dynamics around the sloped CoIn shows a fast rise of the C=CH₂ bond beyond 1.6Å leading to the first spawning event that transfers a significant portion of the population (Figure 8C). Subsequently, after a reduction of the bond length, both the original parent and first child TBF undergo another spawning event. The first child then transfers population to both, the TBF it spawned back to the S₁ as well as the second TBF on the ground state. The two TBFs on the excited state as well as the TBFs on the ground state exchange population between each other until after around 25 fs the original parent TBF spawns again to the ground state leading to the final deactivation of the population. Here, the nature of the spawning algorithm is clearly illustrated as each newly spawned TBF follows its own path on the configuration space driven by the force of the respective electronic state.

The CASPT2 dynamics are all governed by a significantly smaller number of spawning events. They all show an initial elongation of the C=CH₂ bond beyond 1.5Å that also leads to the first spawning event; in contrast to CASSCF the spawn occurs however after the maximum elongation is reached. The population transfer to the spawned child differs greatly between the different CASPT2 flavors being more significant in MS and RMS-CASPT2, as also seen in figure 5. However, while in RMS-CASPT2, the child immediately transfers back population to the parent before independently evolving on the S₀, in MS-CASPT2 the initial parent spawns another child TBF after a few femtoseconds and the two TBFs on the ground state evolve in each others vicinity in configuration space and exchange population with each other. Overall, the evolution in configuration space appears very similar between all CASPT2 flavors despite their large differences in population decay. While the CASSCF bundle evolves similarly along the C=CH₂ bond, the ground state TBFs undergo a much stronger torsion than the CASPT2 ground state. It is interesting to notice that despite the CoIn with higher
torsion angles being predicted as having peaked topology, the evolution of the torsion angles
between the branches of the two initial conditions show only minor differences in the torsion
but seem to mainly differ in the C=CH₂ bond elongation.

<figure><img src="image_9.png" /><figcaption>Figure 8: Time evolution of the C=CH₂ bond length (A and C) and CH₂ torsion angle (B and D) for the same two initial conditions discussed in Figure 6. The evolution of TBFs on the excited state is shown by darker circles, and on the ground state by lighter color with black outlines. The size of the circles is proportional to the population of each TBF of the branch at that point in time.</figcaption></figure>

We here showed how CASSCF and CASPT2 dynamics can differ quite drastically in
the ultrafast time scale. Small variations in how the electronic energies are computed can
also induce relevant differences that might affect the computation of spectroscopic signals,
quantum yields, or kinetic profiles. Our details analysis, easily reproducible thanks to the
user-friendly module available in Py spawn, explains at the same time the reason why these
types of studies are required and why an interface like Py Spawn/OpenMolcas is needed and
essential to have an accurate and balanced description of nuclear dynamics and electronic
structure to disclose the ultrafast dynamics of molecules unambiguously. We hope that this
investigation together with the free availability of our interface will push even more research
efforts in this direction.

Conclusions

In this work, we introduced a brand-new interface between Py spawn, software for running
nonadiabatic molecular dynamics with AIMS algorithm, and OpenMolcas, a package for
running multireference electronic structure calculations. This interface allows the user to
run AIMS/CASSCF and AIMS/CASPT2 dynamics at different flavors: MS-CASPT2, XMS-
CASPT2, XDW-CASPT, and RMS-CASPT2. It exploits the best of the two packages: a
Python, modular implementation of the AIMS algorithm and extensive analysis capability in
Py Spawn, and advanced implementation of electronic structure calculations in OpenMolcas,
including the availability of analytical gradients at the CASPT2 level, crucial for an affordable
dynamics at such an accurate, but costly, level of theory. Though many assume that a swarm
of AIMS calculations are more expensive than a comparable swarm of trajectory surface
hopping calculations, in practice the cost is arguably similar. This is achievable for two
reasons. First, trajectory surface hopping simulations exist discretely in a single electronic
state, stochastically determined at any point in time. Therefore even for a single initial
position and momentum, multiple surface hopping trajectories must be run with different
random seeds to achieve convergence. In contrast, AIMS trajectories have a continuously
variable population that is deterministic, eliminating the need for such over-sampling. Second,
though formally the cost of AIMS simulations scale quadratically with the number of trajectory
basis functions, in practice the vast majority of centroid calculations may be avoided by
thresholding the overlaps of nuclear basis functions. The result is that AIMS simulations scale nearly linearly and statistics converge more rapidly with the number of trajectories than trajectory surface hopping simulations. Depending on the specific physical situation and approximations used, one or the other may be more expensive.

We applied our interface to fulvene, a challenging and interesting system. After being excited to the first singlet electronic state, fulvene can decay to the ground state through two different channels, promoted by two distinct crossings between the $S_0$ and $S_1$ surfaces. The perturbational correction to the energy, changes importantly the potential energy surfaces, closing one of the deactivation channels and inducing a slower relaxation with less degree of the characteristic reflection from $S_0$ to $S_1$. Additionally, the four CASPT2 flavors do not induce the same dynamics, but slight differences can be observed according to the different implementations. XMS-CASPT2 and RMS-CASPT2 are the two flavors that perform the best in nonadiabatic dynamics with the energies calculated on-the-fly, ensuring smooth potentials and energy conservation. The applications and comparison of CASPT2 flavors' performances in nonadiabatic dynamics is an area yet to be fully explored, and in this work, we gave an idea of its importance by pointing out the dependence of the simulated behaviors on the level of electronic structure employed. Our interface represents an important tool offered to the community to further explore such effects and to generally combine accurate nuclear dynamics with the latest developments in OpenMolcas. Last but not least, the implementation is fully open-source, and it can be downloaded with PySpawn on its Git page at https://github.com/blevine37/pySpawn17.

## Acknowledgments

LMI gratefully acknowledges support from the ANR Q-DeLight project, Grant No. ANR-20-CE29-0014 of the French Agence Nationale de la Recherche. AM and BGL gratefully acknowledge support from the National Science Foundation under grant CHE-1954519 and
from the Institute for Advanced Computational Science. This work used Expanse GPU at the
San Diego Supercomputer Center (SDSC) through allocation CHE140101 from the Advanced
Cyberinfrastructure Coordination Ecosystem: Services and Support (ACCESS) program,
which is supported by National Science Foundation grants #2138259, #2138286, #2138307,
#2137603, and #2138296.
References

(1) Mai, S.; González, L. Molecular photochemistry: recent developments in theory. *Angew. Chem. Int. Ed.* **2020**, *59*, 16832–16846.

(2) Janoš, J.; Slavíček, P. What Controls the Quality of Photodynamical Simulations? *Electronic Structure Versus Nonadiabatic Algorithm. J. Chem. Theory Comput.* **2023**, *19*, 8273–8284.

(3) Curchod, B. F.; Martínez, T. J. Ab initio nonadiabatic quantum molecular dynamics. *Chem. Rev.* **2018**, *118*, 3305–3336.

(4) Nelson, T. R.; White, A. J.; Bjorgaard, J. A.; Sifain, A. E.; Zhang, Y.; Nebgen, B.; Fernandez-Alberti, S.; Mozyrsky, D.; Roitberg, A. E.; Tretiak, S. Non-adiabatic excited-state molecular dynamics: Theory and applications for modeling photophysics in extended molecular materials. *Chem. Rev.* **2020**, *120*, 2215–2287.

(5) Worth, G. A.; Burghardt, I. Full quantum mechanical molecular dynamics using Gaussian wavepackets. *Chem. Phys. Lett.* **2003**, *368*, 502–508.

(6) Worth, G. A.; Robb, M. A.; Burghardt, I. A novel algorithm for non-adiabatic direct dynamics using variational Gaussian wavepackets. *Farad. Disc.* **2004**, *127*, 307–323.

(7) Shalashilin, D. V. Quantum mechanics with the basis set guided by Ehrenfest trajectories: Theory and application to spin-boson model. *J. Chem. Phys.* **2009**, *130*.

(8) Shalashilin, D. V. Nonadiabatic dynamics with the help of multiconfigurational Ehrenfest method: Improved theory and fully quantum 24D simulation of pyrazine. *J. Chem. Phys.* **2010**, *132*.

(9) Martinez, T. J.; Levine, R. First-principles molecular dynamics on multiple electronic states: a case study of NaI. *J. Chem. Phys.* **1996**, *105*, 6334–6341.
(10) Martínez, T. J. Ab initio molecular dynamics around a conical intersection: Li (2p) + H2. *Chem. Phys. Lett.* **1997**, *272*, 139–147.

(11) Martinez, T. J.; Ben-Nun, M.; Levine, R. D. Multi-electronic-state molecular dynamics: A wave function approach with applications. *J. Phys. Chem.* **1996**, *100*, 7884–7895.

(12) Martínez, T. J.; Ben-Nun, M.; Levine, R. D. Molecular collision dynamics on several electronic states. *J. Phys. Chem. A* **1997**, *101*, 6389–6402.

(13) Ben-Nun, M.; Molnar, F.; Lu, H.; Phillips, J.; Martínez, T.; Schulten, K. Quantum dynamics of the femtosecond photoisomerization of retinal in bacteriorhodopsin. *Farad. Disc.* **1998**, *110*, 447–462.

(14) Yu, J. K.; Bannwarth, C.; Hohenstein, E. G.; Martinez, T. J. Ab Initio Nonadiabatic Molecular Dynamics with Hole-Hole Tamm-Dancoff Approximated Density Functional Theory. *J. Chem. Theory Comput.* **2020**, *16*, 5499–5511.

(15) Thongyod, W.; Buranachai, C.; Pengpan, T.; Punwong, C. Fluorescence quenching by photoinduced electron transfer between 7-methoxycoumarin and guanine base facilitated by hydrogen bonds: an in silico study. *Phys. Chem. Chem. Phys.* **2019**, *21*, 16258–16269.

(16) Hudock, H. R.; Levine, B. G.; Thompson, A. L.; Satzger, H.; Townsend, D.; Gador, N.; Ullrich, S.; Stolow, A.; Martinez, T. J. Ab initio molecular dynamics and time-resolved photoelectron spectroscopy of electronically excited uracil and thymine. *J. Phys. Chem. A* **2007**, *111*, 8500–8508.

(17) Wang, K.; McKoy, V.; Hockett, P.; Schuurman, M. S. Time-Resolved Photoelectron Spectra of CS2: Dynamics at Conical Intersections. *Phys. Rev. Lett.* **2014**, *112*.

(18) Glover, W. J.; Mori, T.; Schuurman, M. S.; Boguslavskiy, A. E.; Schalk, O.; Stolow, A.; Martinez, T. J. Excited state non-adiabatic dynamics of the smallest polyene, trans 1,3-butadiene. II. Ab initio multiple spawning simulations. *J. Chem. Phys.* **2018**, *148*.
(19) Schuurman, M. S.; Stolow, A. Dynamics at conical intersections. *Annu. Rev. Phys. Chem.* **2018**, *69*, 427–450.

(20) Wolf, T. J. A.; Sanchez, D. M.; Yang, J.; Parrish, R. M.; Nunes, J. P. F.; Centurion, M.; Coffee, R.; Cryan, J. P.; Guehr, M.; Hegazy, K. et al. The photochemical ring-opening of 1,3-cyclohexadiene imaged by ultrafast electron diffraction. *Nat. Chem.* **2019**, *11*, 504–509.

(21) Silfies, M. C.; Mehmood, A.; Kowzan, G.; Hohenstein, E. G.; Levine, B. G.; Allison, T. K. Ultrafast internal conversion and photochromism in gas-phase salicylideneaniline. *J. Chem. Phys.* **2023**, *159*.

(22) Curchod, B. F.; Rauer, C.; Marquetand, P.; González, L.; Martínez, T. J. Communication: GAIMS—Generalized Ab Initio Multiple Spawning for both internal conversion and intersystem crossing processes. *J. Chem. Phys.* **2016**, *144*.

(23) Fedorov, D. A.; Pruitt, S. R.; Keipert, K.; Gordon, M. S.; Varganov, S. A. Ab Initio Multiple Spawning Method for Intersystem Crossing Dynamics: Spin-Forbidden Transitions between 3B1 and 1A1 States of GeH2. *J. Phys. Chem. A* **2016**, *120*, 2911–2919.

(24) Kim, J.; Tao, H.; Martinez, T. J.; Bucksbaum, P. Ab initio multiple spawning on laser-dressed states: a study of 1,3-cyclohexadiene photoisomerization via light-induced conical intersections. *J. Phys. B* **2015**, *48*.

(25) Mignolet, B.; Curchod, B. F.; Martínez, T. J. Communication: XFAIMS—eXternal Field Ab Initio Multiple Spawning for electron-nuclear dynamics triggered by short laser pulses. *J. Chem. Phys.* **2016**, *145*.

(26) Rana, B.; Hohenstein, E. G.; Martinez, T. J. Simulating the Excited-State Dynamics of Polaritons with Ab Initio Multiple Spawning. *J. Phys. Chem. A* **2023**, *128*, 139–151.
(27) Ben-Nun, M.; Martínez, T. J. A multiple spawning approach to tunneling dynamics. *J. Chem. Phys.* **2000**, *112*, 6113–6121.

(28) Virshup, A. M.; Punwong, C.; Pogorelov, T. V.; Lindquist, B. A.; Ko, C.; Martinez, T. J. Photodynamics in complex environments: ab initio multiple spawning quantum mechanical/molecular mechanical dynamics. *J. Phys. Chem. B* **2009**, *113*, 3280–3291.

(29) Song, C. State averaged CASSCF in AMOEBA polarizable water model for simulating nonadiabatic molecular dynamics with nonequilibrium solvation effects. *J. Chem. Phys.* **2023**, *158*.

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

(31) Mignolet, B.; Curchod, B. F. A walk through the approximations of ab initio multiple spawning. *J. Chem. Phys.* **2018**, *148*.

(32) Levine, B. G.; Coe, J. D.; Virshup, A. M.; Martinez, T. J. Implementation of ab initio multiple spawning in the Molpro quantum chemistry package. *Chem. Phys.* **2008**, *347*, 3–16.

(33) Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems. *J. Chem. Theory Comput.* **2013**, *9*, 4959–4972.

(34) Barbatti, M.; Ruckenbauer, M.; Plasser, F.; Pittner, J.; Granucci, G.; Persico, M.; Lischka, H. Newton-X: a surface-hopping program for nonadiabatic molecular dynamics. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2014**, *4*, 26–33.

(35) Akimov, A. V. Libra: An Open-Source “Methodology Discovery” Library for Quantum and Classical Dynamics Simulations. *J. Comput. Chem.* **2016**, *37*, 1626–1649.

(36) 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.

(37) Sifain, A. E.; Bjorgaard, J. A.; Nelson, T. R.; Nebgen, B. T.; White, A. J.; Gifford, B. J.; Gao, D. W.; Prezhdo, O. V.; Fernandez-Alberti, S.; Roitberg, A. E. et al. Photoexcited Nonadiabatic Dynamics of Solvated Push-Pull π-Conjugated Oligomers with the NEXMD Software. *J. Chem. Theory Comput.* **2018**, *14*, 3955–3966.

(38) Mai, S.; Marquetand, P.; Gonzalez, L. Nonadiabatic dynamics: The SHARC approach. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2018**, *8*.

(39) Avagliano, D.; Bonfanti, M.; Nenov, A.; Garavelli, M. Automatized protocol and interface to simulate QM/MM time-resolved transient absorption at TD-DFT level with COBRAMM. *J. Comput. Chem.* **2022**, *43*, 1641–1655.

(40) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O. Second-order perturbation theory with a complete active space self-consistent field reference function. *J. Chem. Phys.* **1992**, *96*, 1218–1226.

(41) Pulay, P. A perspective on the CASPT2 method. *Int. J. Quantum Chem.* **2011**, *111*, 3273–3279.

(42) Malmqvist, P. Å.; Rendell, A.; Roos, B. O. The restricted active space self-consistent-field method, implemented with a split graph unitary group approach. *J. Phys. Chem.* **1990**, *94*, 5477–5482.

(43) Borrego-Varillas, R.; Nenov, A.; Kabaciński, P.; Conti, I.; Ganzer, L.; Oriana, A.; Jaiswal, V. K.; Delfino, I.; Weingart, O.; Manzoni, C. et al. Tracking excited state decay mechanisms of pyrimidine nucleosides in real time. *Nat. Comm.* **2021**, *12*, 7285.

(44) Park, J. W.; Shiozaki, T. On-the-fly CASPT2 surface-hopping dynamics. *J. Chem. Theory Comput.* **2017**, *13*, 3676–3683.
(45) Polyak, I.; Hutton, L.; Crespo-Otero, R.; Barbatti, M.; Knowles, P. J. Ultrafast photoin-
duced dynamics of 1, 3-cyclohexadiene using XMS-CASPT2 surface hopping. *J. Chem.
Theory Comput.* **2019**, *15*, 3929–3940.

(46) Coe, J. D.; Levine, B. G.; Martínez, T. J. Ab initio molecular dynamics of excited-state intramolecular proton transfer using multireference perturbation theory. *J. Phys. Chem. A* **2007**, *111*, 11302–11310.

(47) Barneschi, L.; Kaliakin, D.; Huix-Rotllant, M.; Ferré, N.; Filatov, M.; Olivucci, M. Assessment of the Electron Correlation Treatment on the Quantum-Classical Dynamics of Retinal Protonated Schiff Base Models: XMS-CASPT2, RMS-CASPT2, and REKS Methods. *J. Chem. Theory Comput.* **2023**, *19*, 8189–8200.

(48) Heindl, M.; González, L. A XMS-CASPT2 non-adiabatic dynamics study on pyrrole.
*Comput. Theor. Chem.* **2019**, *1155*, 38–46.

(49) Nishimoto, Y. Analytic gradients for restricted active space second-order perturbation theory (RASPT2). *J. Chem. Phys.* **2021**, *154*.

(50) Nishimoto, Y.; Battaglia, S.; Lindh, R. Analytic first-order derivatives of (X) MS, XDW, and RMS variants of the CASPT2 and RASPT2 methods. *J. Chem. Theory Comput.* **2022**, *18*, 4269–4281.

(51) Celani, P.; Werner, H.-J. Analytical energy gradients for internally contracted second-order multireference perturbation theory. *J. Chem. Phys.* **2003**, *119*, 5044–5057.

(52) MacLeod, M. K.; Shiozaki, T. Communication: Automatic code generation enables nuclear gradient computations for fully internally contracted multireference theory. *J. Chem. Phys.* **2015**, *142*.

(53) Song, C.; Neaton, J. B.; Martínez, T. J. Reduced scaling formulation of CASPT2 analytical gradients using the supporting subspace method. *J. Chem. Phys.* **2021**, *154*.
(54) Shiozaki, T.; Győrffy, W.; Celani, P.; Werner, H.-J. Communication: Extended multi-state complete active space second-order perturbation theory: Energy and nuclear gradients. *J. Chem. Phys.* **2011**, *135*.

(55) Battaglia, S.; Lindh, R. Extended dynamically weighted CASPT2: The best of two worlds. *J. Chem. Theory Comput.* **2020**, *16*, 1555–1567.

(56) Battaglia, S.; Lindh, R. On the role of symmetry in XDW-CASPT2. *J. Chem. Phys.* **2021**, *154*.

(57) Fdez. Galván, I.; Vacher, M.; Alavi, A.; Angeli, C.; Aquilante, F.; Autschbach, J.; Bao, J. J.; Bokarev, S. I.; Bogdanov, N. A.; Carlson, R. K. et al. OpenMolcas: From source code to insight. *J. Chem. Theory Comput.* **2019**, *15*, 5925–5964.

(58) Li Manni, G.; Fdez. Galván, I.; Alavi, A.; Aleotti, F.; Aquilante, F.; Autschbach, J.; Avagliano, D.; Baiardi, A.; Bao, J. J.; Battaglia, S. et al. The OpenMolcas Web: A community-driven approach to advancing computational chemistry. *J. Chem. Theory Comput.* **2023**, *19*, 6933–6991.

(59) 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.

(60) Meek, G. A.; Levine, B. G. The best of both Reps—Diabatized Gaussians on adiabatic surfaces. *J. Chem. Phys.* **2016**, *145*, 184103.

(61) Finley, J.; Malmqvist, P.-Å.; Roos, B. O.; Serrano-Andrés, L. The multi-state CASPT2 method. *Chem. Phys. Lett.* **1998**, *288*, 299–306.

(62) Battaglia, S.; Lindh, R. On the role of symmetry in XDW-CASPT2. *J. Chem. Phys.* **2021**, *154*.
(63) Schreiber, M.; Silva-Junior, M. R.; Sauer, S.; Thiel, W. Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3. *J. Chem. Phys.* **2008**, *128*.

(64) Granovsky, A. A. Extended multi-configuration quasi-degenerate perturbation theory: The new approach to multi-state multi-reference perturbation theory. *J. Chem. Phys.* **2011**, *134*.

(65) Zobel, J. P.; Nogueira, J. J.; González, L. The IPEA dilemma in CASPT2. *Chem. Sci.* **2017**, *8*, 1482–1499.

(66) Sauri, V.; Serrano-Andrés, L.; Shahi, A. R. M.; Gagliardi, L.; Vancoillie, S.; Pierloot, K. Multiconfigurational second-order perturbation theory restricted active space (RASPT2) method for electronic excited states: A benchmark study. *J. Chem. Theory Comput.* **2011**, *7*, 153–168.

(67) 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.

(68) 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.

(69) Alfalah, S.; Belz, S.; Deeb, O.; Leibscher, M.; Manz, J.; Zilberg, S. Photoinduced quantum dynamics of ortho-and para-fulvene: Hindered photoisomerization due to mode selective fast radiationless decay via a conical intersection. *J. Chem. Phys.* **2009**, *130*.

(70) 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*.
(71) Weight, B. M.; Mandal, A.; Huo, P. Ab initio symmetric quasi-classical approach to investigate molecular Tully models. *J. Chem. Phys.* **2021**, *155*.

(72) Ibele, L. M.; Lassmann, Y.; Martínez, T. J.; Curchod, B. F. Comparing (stochastic-selection) ab initio multiple spawning with trajectory surface hopping for the photodynamics of cyclopropanone, fulvene, and dithiane. *J. Chem. Phys.* **2021**, *154*.

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

(74) 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.

(75) Weight, B.; Mandal, A.; Hu, D.; Huo, P. Ab initio Spin-Mapping Non-adiabatic Dynamics Simulations of Photochemistry. 2024; ChemRxiv; doi:10.26434/chemrxiv-2024-4hzlj.
TOC Graphic

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