Comparing (stochastic-selection) ab initio multiple spawning with trajectory surface hopping for the photodynamics of cyclopropanone, fulvene, and dithiane

Ab Initio Multiple Spawning (AIMS) simulates the excited-state dynamics of molecular systems by representing nuclear wavepackets in a basis of coupled traveling Gaussian functions, called trajectory basis functions (TBFs). New TBFs are spawned when nuclear wavepackets enter regions of strong nonadiabaticity, permitting the description of non-Born–Oppenheimer processes. The spawning algorithm is simul-taneously the blessing and the curse of the AIMS method: it allows for an accurate description of the transfer of nuclear amplitude between different electronic states, but it also dramatically increases the computational cost of the AIMS dynamics as all TBFs are coupled. Recently, a strategy coined stochastic-selection AIMS (SSAIMS) was devised to limit the ever-growing number of TBFs and tested on simple molecules. In this work, we use the photodynamics of three different molecules—cyclopropanone, fulvene, and 1,2-dithiane—to investigate (i) the potential of SSAIMS to reproduce reference AIMS results for challenging nonadiabatic dynamics, (ii) the compromise achieved by SSAIMS in obtain-ing accurate results while using the smallest average number of TBFs as possible, and (iii) the performance of SSAIMS in comparison to the mixed quantum/classical method trajectory surface hopping (TSH)—both in terms of its accuracy and computational cost. We show that SSAIMS can accurately reproduce the AIMS results for the three molecules considered at a much cheaper computational cost, often close to that of TSH. We deduce from these tests that an overlap-based criterion for the stochastic-selection process leads to the best agreement with the reference AIMS dynamics for the smallest average number of TBFs.


I. INTRODUCTION
Simulating the dynamics that molecules undergo after light absorption poses a complete challenge for theoretical chemistry as this implies moving beyond the celebrated Born-Oppenheimer approximation. 1,2 Following photoexcitation, molecules are likely to access regions of configuration space where nuclear motion can trigger changes in electronic eigenstates-the so-called nonadiabatic effects-causing a breakdown of the Born-Oppenheimer approximation. As a result, a solution of the coupled electronnuclear time-dependent Schrödinger equation is required to investigate the photophysical and photochemical processes of molecules. However, this is often hampered by the dimensionality of the problem when looking at molecules-introducing approximations becomes inevitable.
A typical starting point for developing nonadiabatic molecular dynamics techniques is to express the molecular wavefunction within the Born-Huang representation. This representation introduces the common picture of photochemical processes where nuclear wavefunctions evolve on potential energy surfaces and transfer amplitudes between different electronic states during nonradiative relaxation processes. Highly accurate methods, such as multiconfigurational time-dependent Hartree (MCTDH), [3][4][5][6] express the electronic structure quantities and nuclear wavefunctions on a grid and subsequently allow for a numerically exact solution of the nuclear time-dependent Schrödinger equation for a few tens of nuclear degrees of freedom. MCTDH is widely viewed as the goldstandard for nonadiabatic dynamics methods. However, only a subset of the nuclear coordinates can be considered when simulating the nonadiabatic dynamics of larger molecules due to the computational cost of the technique and the need for precomputed potential energy surfaces (see Refs. 7 and 8 for recent developments on this topic). Thus, one should consider additional approximations to the molecular time-dependent Schrödinger equation to simulate the excited-state dynamics of molecules in their full configurational space.
One of the most famous families for nonadiabatic dynamics is the so-called mixed quantum/classical approach, where nuclei are treated classically and electrons quantum-mechanically. Their simplicity and comparably low cost make mixed quantum/classical methods very appealing, despite their underlying approximations. Techniques such as Ehrenfest dynamics and trajectory surface hopping (TSH) have evolved to become popular choices for nonadiabatic dynamics in many cases. 9,10 By applying an independent trajectory approximation (ITA), TSH mimics the dynamics of nuclear wavepackets by propagating a swarm of totally independent classical trajectories; nonadiabatic effects are included by allowing the trajectories to hop between electronic states stochastically (we will focus in this work only on the fewest-switches flavor of TSH). 9,11-13 While TSH has been successfully employed to describe the photochemistry and photophysics of a great variety of molecular systems in the past, its formalism limits the possible improvements to some of its deficiencies, such as the overcoherence problem. More details on TSH will be given below. More recent developments of mixed quantum/classical strategies proposed to partly alleviate the independent trajectory approximation by including some form of couplings between trajectories. 14,15 Another family of methods for nonadiabatic dynamics proposes to express the nuclear wavefunctions in a linear combination of coupled trajectory basis functions (TBFs), which are nothing but multidimensional Gaussian functions traveling in configuration space. A broad variety of TBF-based methods have been developed over the past decades, such as variational multiconfigurational Gaussian (vMCG), [16][17][18] multiconfiguration Ehrenfest (MCE), [19][20][21] and full and ab initio multiple spawning (FMS and AIMS). [22][23][24] The differences between these strategies reside in the way the TBFs are propagated and the level of approximations applied to the description of their mutual couplings. vMCG employs a quantum propagation of the TBFs based on the Dirac-Frenkel variational principle, while MCE, FMS, and AIMS propagate the TBFs classically. MCE uses a mean-field potential energy for the propagation of the TBFs. FMS and AIMS propagate the TBFs adiabatically on the potential energy surfaces but allow for an expansion of the basis set to describe nonadiabatic transitions through the so-called spawns (see Sec. II A for more details). In stark contrast with TSH, all these methods are formally exact and can be derived from first-principles. Their accuracy is only limited by the number of TBFs used to expand the nuclear wavefunctions and the description of the couplings between the TBFs.
In theory, the number of electronic structure calculations required per AIMS time step scales quadratically with the number of TBFs due to their mutual coupling. In most implementations of the AIMS algorithm, however, this scaling is lowered by neglecting coupling elements between TBFs with vanishing nuclear overlap. 25 In addition, the spawning algorithm continually increases the size of the basis when encountering regions of nonadiabaticity. Hence, in cases of repeated nonadiabatic transitions, especially if more than two electronic states are considered in the dynamics, the number of TBFs and thus the overall cost of the dynamics can increase dramatically. As the spawning algorithm is only designed to create new TBFs, it is not unusual that a large number of TBFs are carried throughout the dynamics with the nuclear amplitude only distributed among a few of them-the other TBFs do not contribute to the description of the nuclear wavepackets anymore or are not coupled to the rest of the swarm. Recently, an algorithm was proposed to reduce the number of TBFs propagated in AIMS by exploiting the fact that groups of TBFs can become uncoupled from each other in the course of the dynamics, allowing to devise selection rules for keeping TBFs alive (see Sec. II A for details). This approach, coined stochastic-selection ab initio spawning (SSAIMS), 26 offers a way to reduce the cost of an AIMS simulation with only a minimal loss of accuracy in the description of the nonadiabatic processes. Proof-of-principle tests of SSAIMS have revealed its potential as a cheaper alternative to AIMS. In this work, we propose to put SSAIMS under a stringent test by simulating the challenging excitedstate dynamics of three molecules-cyclopropanone, fulvene, and 1,2-dithiane-and compare its accuracy and computational cost to both AIMS and TSH.
To motivate the adoption of cyclopropanone, fulvene, and 1,2dithiane as a testbed for SSAIMS, we provide here a brief description of their photodynamics, illustrated by a schematic depiction of their nonradiative decay (Fig. 1). The dynamics of cyclopropanone has been subject to several studies in the past. [27][28][29][30] Photoexcitation to the first excited singlet state (S 1 ) triggers a ring-opening reaction, mediated by carbon-carbon bond breaking followed by subsequent dissociation into ethylene and CO. The S 0 and S 1 potential energy curves along an interpolation in internal coordinates between the ring-closed Franck-Condon and ring-opened minimum energy conical intersection are depicted schematically in Fig. 1(a). The two electronic states come close in energy during the ring-opening process and separate after the nonadiabatic transition to the ground state. Thus, for cyclopropanone, one expects a relatively straightforward transfer of the nuclear wavepacket from S 1 to S 0 without a significant contribution of population back transfer to S 1 . However, the nonradiative decay is governed by strong geometrical deformations of the molecule-ring opening and subsequent dissociation-which can challenge the employed electronic structure methods. During an AIMS simulation, it is sufficient that one of the TBFs suffers from an electronic structure instability to stop the propagation of all the TBFs, even if these are only weakly coupled. SSAIMS could reduce the probability of such a dramatic issue by ensuring that only strongly coupled TBFs are propagated together. Following photoexcitation, fulvene can undergo a reflection process at a strongly sloped conical intersection between S 1 and S 0 31-34 and has subsequently been compared to a molecular version of Tully's third model. 35 Figure 1(b) shows a one-dimensional cut of the S 0 and S 1 potential energy surfaces in the vicinity of the sloped conical intersection. Coming from the Franck-Condon region in S 1 , the nuclear wavepacket crosses the sloped conical intersection and, soon after being reflected on S 0 , splits into a part decaying to the S 0 Franck-Condon region and a part transferring back to the S 1 state after meeting the same conical intersection. A proper description of this process must account for the interaction between the populations on the different excited states. As a result, the excited-state dynamics of fulvene appeared to be highly sensitive to the nonadiabatic dynamics methods employed 35 and therefore provides an interesting test for SSAIMS. Finally, 1,2-dithiane shows a sulfur-sulfur ring-opening process upon photoexcitation in its first excited electronic state, followed by a complex dynamics where the opened ring recloses within 300 fs. 36 By looking at the three lowest singlet states along a linear interpolation in internal coordinates between the Franck-Condon region, the S 1 minimum geometry, and the S 1 /S 0 minimum energy conical intersection [ Fig. 1(c)], it becomes apparent that the three singlet states become (and remain) nearly degenerate soon after the ring opening. An accurate theoretical description of the dynamics poses a challenge because of the complex ringopening and ring-closing process. Furthermore, the near degeneracy of the three lowest singlet states can induce repeated population transfer-yet another challenge for nonadiabatic dynamics methods.
In the following, we present the different methods for nonadiabatic dynamics that will be compared in this work and propose an extensive discussion of our computational details (Sec. II). We then assess the performance of AIMS, TSH, and SSAIMS for the three different molecules discussed above in Sec. III, focusing on the population trace over time as well as the computational effort of each method. Our main finding is that Overlap SSAIMS (OSSAIMS), one flavor of SSAIMS, not only reproduces the results obtained with AIMS for all three molecules, but it does so at a similar computational cost than TSH.

A. Theoretical methods
In the following, we offer a brief introduction and overview of the key methods employed in this work-TSH, AIMS, and SSAIMS.
The reader is referred to more general reviews on AIMS, 37,38 TSH, 13,39 and nonadiabatic dynamics in general 40 for additional details.

Trajectory surface hopping
Due to its mixed quantum/classical nature, TSH constitutes an appealing approach to nonadiabatic dynamics. Its underlying independent trajectory approximation (ITA) allows for the use of classical, fully independent trajectories to describe the nuclear wavepacket dynamics. This, in turn, leads to a dramatic reduction of complexity in describing the dynamics of nuclear wavepackets and, as a result, a lowering of the computational cost. In a typical TSH dynamics, one starts by choosing a set of initial conditions to represent the nuclear wavepacket at time t = 0. Each independent trajectory is then propagated using a classical nuclear force obtained from the potential energy surface driving the dynamics. The probability of hopping to a different electronic state is based on the strength of the corresponding nonadiabatic couplings. Based on the calculated hopping probabilities, a stochastic process indicates, after each time step, whether a surface hop should occur. The most commonly used version of TSH is the fewest switches algorithm, proposed by Tully in 1990. 11 In the following, "TSH" is used synonymously to "fewest switches surface hopping." While computationally appealing, the ITA limits the accuracy of TSH as it hampers a proper description of the branching of a nuclear wavepacket following a nonadiabatic transition-when the different forces coming from the coupled electronic states split the wavepacket in configuration space. This issue is sometimes called the decoherence problem of TSH, [41][42][43][44][45][46][47][48] and a variety of ad hoc corrections have been proposed to alleviate this issue. 39,45,46,[49][50][51][52] 2. Ab initio multiple spawning AIMS is an approach derived from the in principle exact framework of full multiple spawning (FMS) to perform on-the-fly nonadiabatic quantum molecular dynamics. [22][23][24] In FMS, a basis of moving multidimensional Gaussian functions, called TBFs, is used to expand the nuclear amplitudes in the Born-Huang expansion. Each TBF is associated with a complex amplitude, and the swarm of TBFs can therefore be pictured as a form of moving grid that supports the propagation of the nuclear wavefunctions. The TBFs evolve classically, following a given adiabatic potential energy surface. When a TBF encounters a region where the nonadiabatic coupling between its running state and another electronic state exceeds a certain predefined threshold, it can spawn a new child TBF on the coupled state. The two TBFs then evolve in a fully coupled manner, which allows for a smooth transfer of amplitude between the two electronic states. Importantly, all the TBFs present during the dynamics are coupled through the time-dependent Schrödinger equation so that both inter-and intrastate couplings are considered, allowing for the exchange of nuclear amplitudes between TBFs. In the limit of a sufficiently large number of TBFs, FMS would become formally exact [see Ref. 53 for a numerical demonstration]. However, the exact treatment of the coupling terms between TBFs would require the knowledge of the potential energy surfaces and nonadiabatic coupling terms over the entire nuclear configuration space, i.e., a precomputation of all these quantities. Therefore, an on-thefly propagation of the FMS equations for molecules in their full dimensionality requires two central approximations, leading to the method called AIMS. The saddle-point approximation (SPA) proposes to take advantage of the spatial localization of the TBFs to approximate the electronic energies and the nonadiabatic couplings using a Taylor series (to zeroth order) around the centroid position between the two coupled TBFs. This strategy dramatically simplifies the evaluation of the intra-and interstate coupling terms between TBFs. In the independent first generation approximation (IFGA), the initial TBFs that represent the nuclear wavefunction at time t = 0 are considered uncoupled during the dynamics. These two controlled approximations constitute the framework of the AIMS method, and their effect on the dynamics has recently been tested. 53,54

Stochastic-selection ab initio multiple spawning
While the SPA and IFGA introduced in the previous paragraph make AIMS computationally tractable for molecules, the cost of this technique still formally scales quadratically with the number of TBFs (N TBF ), N TBF × (N TBF + 1)/2 to be precise. This scaling can hamper the use of AIMS for molecules that encounter repeated crossings between surfaces, as this would lead to a dramatic increase in the number of TBFs. The Achilles' heel of AIMS in its original formulation is that spawns increase the number of TBFs, but no systematic death process exists to control this ever-growing population. A new algorithm has been proposed recently to address this issue and leverages the naturally occurring decoupling between groups of TBFs. 26 By stochastically selecting subgroups of coupled TBFs, this so-called stochastic-selection ab initio multiple spawning (SSAIMS) offers an efficient way of reducing the cost of typical AIMS dynamics within a controlled framework of approximations. An SSAIMS run is very similar to an AIMS one (see Fig. 2 for a schematic representation): An initial (parent) TBF evolves on its respective electronic state and spawns new TBFs in regions of significant nonadiabatic coupling, which in turn can spawn other TBFs at will (new lines and TBFs in Fig. 2). The essential difference is that in SSAIMS, one monitors the coupling between all TBFs at every time step and tries to detect when TBFs (or groups of TBFs) become uncoupled. This uncoupling occurs when the coupling drops below a certain threshold ϵ. If uncoupled (groups of) TBFs are detected, the SSAIMS algorithm takes effect: one computes the population of all uncoupled groups of TBFs, a random number is generated, and one selects one of the groups of TBFs to continue the simulation based on a Monte Carlo procedure ("Stochastic selection" in Fig. 2). The unselected TBFs are stopped at this point of the dynamics, and the population of The coupling between TBFs is constantly monitored during the dynamics, and when TBFs become uncoupled, a stochastic-selection process is triggered. Only one TBF, or a group of coupled TBFs, will survive the stochastic-selection process-In the present scheme, the parent TBF becomes uncoupled from a group of two child TBFs and is therefore discarded.
the remaining TBFs is renormalized. By repeating the same initial condition with different seeds for the random number generator, one can show that the SSAIMS result converges toward the AIMS one-a strategy called SSAIMS with repetitions.
In AIMS, the TBFs rapidly spread in phase space due to their classical propagation and become decoupled so that many independent, redundant TBFs may be carried along the simulation. SSAIMS detects those and removes them from the simulation based on a stochastic selection process. Thus, in the limit of a minimum selection threshold and a sufficient number of repetitions for each initial condition, SSAIMS will converge toward AIMS, as demonstrated in Ref. 26 (the interested reader is also referred to this work for a detailed presentation of SSAIMS). The coupling between two TBFs can be defined in different ways, and in the following, we will focus on two quantities: the absolute value of the Hamiltonian matrix elements between the TBFs considered (Energy-SSAIMS, ESSAIMS) and the absolute value of the overlap between the two TBFs (Overlap-SSAIMS, OSSAIMS).
The stochastic nature of SSAIMS makes this method somehow reminiscent of some of the ideas of TSH. However, the stochastic processes at play are different in these two methods. In TSH, the stochastic process determines the nonadiabatic transfer of population per se. In contrast, the nonadiabatic process in SSAIMS is described exactly as in AIMS, as the coupling between the TBFs is maximal in this part of the dynamics, and the stochastic-selection process will have no effect. It is only after the branching of the TBFs that uncoupling may start to appear, and stochastic-selection processes can occur. As mentioned above, TSH relies on ad hoc corrections to approximate decoherence effects. In contrast, AIMS accounts for decoherence intrinsically by virtue of the coupling between TBFs, which evolve based on the nuclear forces determined for their respective electronic state. Importantly, this natural account for decoherence is not compromised in SSAIMS.

Computational cost of AIMS, SSAIMS, and TSH
This work aims not only to compare the performance of SSAIMS in describing complex nonadiabatic dynamics processes with that of AIMS and TSH but also to analyze their respective computational costs.
A comparison of the computational cost in terms of electronic structure calls or wall time is plagued by the algorithmic details related to the implementation of the method (e.g., adaptive time steps, convergence criteria, or propagation algorithms), which hampers a formal comparison of the computational burden associated with each technique. To provide a comparison that is as fair as possible, we propose here to focus on the "theoretical number of electronic structure calculations" required at each time step of the dynamics, which unravels the computational effort of each method for the different molecular systems presented. In TSH, this number remains constant over time as only one electronic structure call is necessary for each trajectory at every time step, i.e., it is simply the product of number of initial conditions (N init ) with the number of runs per initial condition (nrun), N TSH ES = nrun × N init . In contrast, the number of electronic structure calls in AIMS and (E/O)SSAIMS will depend on the number of TBFs present in the dynamics at every time step. In the following, we will consider the worst case scenario for AIMS and (E/O)SSAIMS, in which coupling elements between all TBFs are always calculated during the dynamics. In this limit, the theoretical number of electronic structure calculations for AIMS is given by N AIMS needs to further account for the number of runs per initial condition (and the fact that the number of TBFs created for each initial condition might not be the same within different runs), that is, In practice, an AIMS calculation can be made cheaper by monitoring the overlap between TBFs to determine whether it is necessary to compute their intra-or interstate coupling terms, as detailed in Ref. 25. The theoretical number of electronic structure calculations presented in this work is, without a doubt, an idealized representation of the computational effort of TSH and (E/O)SSAIMS. Still, it offers a formal and straightforward way of comparing the cost of these different methods and unraveling the complexity engendered by the couplings between TBFs in AIMS-based techniques.

Electronic structure
Electronic energies, nuclear gradients of the energies, and nonadiabatic couplings were computed for all molecules studied in this work with state-averaged complete active space selfconsistent field (SA-CASSCF) 55,56 and a 6-31G * basis set. 57,58 For cyclopropanone, an (8/7) active space was used, consisting of the C= =O ππ * orbitals, the σσ * pairs of the adjacent C− −C bonds, and one n lone pair orbital of the oxygen atom, and the state averaging procedure was performed for the lowest two singlet states. For fulvene, a state averaging for the lowest two singlet states and a (6/6) active space were employed, including the three pairs of ππ * orbitals. SA-CASSCF calculations for cyclopropanone and fulvene were carried out with the MOLPRO 2012 program package. 59 The electronic structure of 1,2-dithiane was described with a (6/4) active space [see the supplementary material for a comparison with a larger (10/8) active space] that includes the σσ * pair of the S-S bond as well as two n lone pair orbitals on the sulfur atoms and with a state averaging for the lowest three singlet states. The SA-CASSCF calculations for 1,2-dithiane were performed with the MOLPRO 2012 program package 59 (for the TSH dynamics) and TeraChem 60-64 for the AIMS and (E/O)SSAIMS dynamics.

Nuclear dynamics
The AIMS and (E/O)SSAIMS dynamics for cyclopropanone and fulvene were performed with the AIMS/MOLPRO interface. 25 For 1,2-dithiane, the AIMS implementation in TeraChem was used. 65,66 All AIMS and (E/O)SSAIMS dynamics share the very same set of parameters for each molecule. The TBFs were propagated with a time step of 20 atomic time units (atu), reduced to 5 atu in regions of strong nonadiabaticity. The criterion to enter the spawning mode used the norm of the nonadiabatic coupling vectors, with a threshold set to 3.0 a.u. −1 for cyclopropanone, 10.0 a.u. −1 for fulvene, and 20.0 a.u. −1 for 1,2-dithiane. The minimum population required for a TBF to spawn was 0.01 for cyclopropanone and fulvene and 0.05 for 1,2-dithiane. The maximum acceptable overlap between a newly spawned TBF and the existing pool of TBFs was set to 0.6 for fulvene and cyclopropanone and 0.5 for 1,2-dithiane. For cyclopropanone, the TBFs running on the ground electronic state were stopped if they remained uncoupled with any other TBFs for more than 100 atu. The threshold for total (classical) energy conservation was set to 0.019 a.u. for cyclopropanone, 0.01 a.u. for fulvene, and 0.005 a.u. for 1,2-dithiane. The same set of initial conditions were used for both (E/O)SSAIMS and AIMS, but the (E/O)SSAIMS runs for each initial condition were repeated five times (with different initial random seeds) to converge the stochastic processes. Different thresholds for the stochastic selection were used and are detailed below or discussed in the supplementary material.
All TSH simulations were performed with the SHARC 2.0 program. [67][68][69] The TSH trajectories were initiated from the very same set of initial conditions as the AIMS/SSAIMS parent TBFs. Every trajectory was repeated multiple times, typically between five and eight times, with different random seeds to ensure convergence of the stochastic process for the nonadiabatic transitions. The number of repetitions was determined such that the maximum standard error of the S 1 population decay of ESSAIMS and (decoherencecorrected) TSH was in agreement: This criterion was fulfilled with five runs for cyclopropanone, seven for fulvene, and eight for 1,2dithiane (for each initial condition). The integration time step for the nuclear propagation was set to 0.5 fs [to resemble 20 atu used in (E/OSS)AIMS], a local diabatization scheme was used, and the nonadiabatic couplings were constructed from wavefunction overlaps to avoid the explicit computation of the nonadiabatic coupling vectors. 70 Following a surface hop, the nuclear kinetic energy was rescaled along the nuclear momenta.
All TSH simulations were carried out with and without the energy-based decoherence correction (EDC) scheme, 49,71 which accounts for decoherence through a dampening of the electronic populations in TSH. The decoherence parameter C was set to 0.1 a.u., as proposed in Ref. 49. The default implementation of the EDC in SHARC was applied, and the same random seeds were used for both TSH and TSH with a decoherence correction (called dTSH of Chemical Physics ARTICLE scitation.org/journal/jcp in the following). We note that the original implementation of the EDC in SHARC was employed, where the amplitudes of the nonrunning states are damped, instead of the populations, as done, for example, in Newton-X. 49,72,73 All the dTSH trajectories strictly satisfy the internal consistency criterion discussed in Ref. 71, and the standard error of the dTSH population decay was calculated using the electronic populations of the trajectories. All the initial conditions employed in this work as well as all the population traces are provided as the supplementary material of this article.

III. RESULTS AND DISCUSSION
A. Cyclopropanone-An exemplar scenario of coupled TBFs leading to instabilities in AIMS The nonadiabatic dynamics of cyclopropanone has been heavily studied in the past employing different combinations of TSH flavors and electronic structure methods. [27][28][29][30] All these previous studies highlight the fact that the S 1 population decay of cyclopropanone following photoexcitation is rather straightforward-the nuclear wavepacket decays from S 1 to S 0 without any significant back population transfer. The S 1 decay has been characterized by a latency time in S 1 of ∼25 fs, followed by a stepwise decay of the S 1 population. 30 However, the structural evolution of the molecule is intriguing: Upon excitation, one or two of the C− −C bonds adjacent to the carbonyl moiety are broken, leading to a dissociation into CO and ethylene within several hundreds of femtoseconds. [27][28][29] Because of these substantial geometrical rearrangements-a ring-opening process, followed by a full dissociation-cyclopropanone poses an interesting challenge for nonadiabatic dynamics methods, in particular for the underlying electronic structure method.

AIMS vs (E/O)SSAIMS
We start by comparing the S 1 population trace obtained by AIMS for the photodynamics of cyclopropanone with that of (E/O)SSAIMS for different selection thresholds (lower panel of Fig. 3). Due to electronic structure instabilities, AIMS (black line in Fig. 3) can only simulate the first 50 fs of dynamics. The active space employed for SA-CASSCF is not stable when the molecule dissociates on the ground state. While the coupling between TBFs in AIMS permits an adequate description of decoherence processes, it comes with the severe drawback that any instabilities in the propagation of one of the TBFs will prevent the propagation of the entire branch of coupled TBFs. In the particular case of cyclopropanone, a dissociating TBF on the ground state remains coupled, even if only weakly, to other TBFs, and any electronic structure instability following this dissociation makes it impossible to run the AIMS dynamics for longer times.
This situation highlights one of the key advantages of the stochastic selection algorithm: SSAIMS can detect when a TBF (here evolving on the ground state) is only weakly coupled to the remaining swarm of TBFs, and perform a selection process accordingly. Applying ESSAIMS with the smallest possible energy threshold (ϵ = 10 −10 a.u.) allows us to prolong the dynamics up to 75 fs (see the light blue dashed curve in Fig. 3). The agreement between AIMS and ESSAIMS (ϵ = 10 −10 a.u.) in the first 50 fs is excellent, with the ESSAIMS population trace remaining well within the standard error of AIMS. Both methods show that the S 1 population starts decaying after 25 fs, before plateauing at 70% of S 1 population. ESSAIMS (ϵ = 10 −10 a.u.) remains stable enough to describe the beginning of the subsequent S 1 population decay. The average number of TBFs for AIMS peaks at 1.4 TBFs (top panel of Fig. 3) and stays near this value until the dynamics fails. Despite a relatively low selection threshold, the average number of TBFs in ESSAIMS (ϵ = 10 −10 a.u.) is already reduced with a peak under 1.25. This number then drops to 1.0 before rising again for the second decay of the S 1 population. These results show that some TBFs in AIMS are likely to be only very weakly coupled. Using a conservative energy threshold for ESSAIMS allows us to stochastically select some weakly coupled TBFs evolving in S 0 but does not fully resolve the instability issues. Increasing the ESSAIMS threshold permits to enforce a faster uncoupling of the TBFs in the dynamics. This approximation should further help with the remaining instabilities.
With a threshold of ϵ = 10 −5 a.u., ESSAIMS is stable enough to describe the full S 1 population decay (see the blue curve in Out of curiosity, one can finally test how ESSAIMS would behave when the selection threshold is set to a very high value, here ϵ = 1 a.u. (see the light blue thin curve in Fig. 3). In this extreme case, the stochastic selection is triggered immediately after the spawn of a new TBF (as can be observed in the average number of TBFs), hindering the exchange of amplitude between the two coupled electronic states.
Let us now test the performance of OSSAIMS, where the criterion for the coupling between TBFs is given by the overlap between the two TBFs considered. An overlap threshold of ϵ = 0.5, which may at first glance appear rather large, leads to a very close agreement with the population decay observed in AIMS and ESSAIMS (dashed light red curve in Fig. 3). Moreover, OSSAIMS achieves this result by requiring a consistently lower average number of TBFs than ESSAIMS (ϵ = 10 −5 a.u.). In fact, OSSAIMS (ϵ = 0.5) uses on average nearly the same number of TBFs as the extreme ESSAIMS (ϵ = 1.0 a.u.) with the major difference that ESSAIMS (ϵ = 1.0 a.u.) leads to a poor description of the nonadiabatic processes. Reducing the OSSAIMS threshold to ϵ = 0.1 does not improve the S 1 population trace significantly, which remains well within the error bars of ESSAIMS (ϵ = 10 −5 a.u.). Once more, the OSSAIMS average number of TBFs with ϵ = 0.1 is consistently lower than for ESSAIMS (ϵ = 10 −5 a.u.), indicating that OSSAIMS could offer a better compromise between accuracy and efficiency.

Comparison between AIMS, (E/O)SSAIMS, and (d)TSH
Now that we showed how (E/O)SSAIMS makes a multiple spawning simulation of cyclopropanone possible, we wish to compare these results with the ones obtained using the mixed quantum/classical method (d)TSH. TSH with and without a decoherence correction (dTSH and TSH, respectively, shown in the lower part of Fig. 4) describes a very similar S 1 population decay for the first 75 fs of dynamics. The population decay starts after 25 fs and plateaus at around 65% of S 1 population. After 50 fs, the population transfer resumes and a difference starts appearing between TSH and dTSH after 70 fs of dynamics, the former showing a slower population decay than the latter. This deviation between the two methods can be rationalized as follows. After a hop to the ground state, the decoherence correction enforces a quenching of the dTSH electronic population to the ground state. In TSH, the electronic coefficient for S 1 is not dampened, artificially increasing the probabilities of hops back to the S 1 state, leading to an overall slowdown of the S 1 population decay.
The initial population decay observed in (d)TSH at 25 fs agrees with (E/O)SSAIMS and AIMS. However, the population trace plateaus in (d)TSH at a lower value than the spawning methods and is outside of the AIMS standard error. dTSH then predicts a faster decay of the S 1 population than ESSAIMS (ϵ = 10 −5 a.u.). While the population decay of TSH seems to agree well with the ESSAIMS population at later times, this is most likely only an artifact of TSH overcoherence (as detailed above). We note that a similar effect of TSH overcoherence was observed in the photodynamics of ethylene. 35 Let us now focus on the computational cost of the different methods compared above (upper panel of Fig. 4), using the theoretical number of electronic structure calls per time step, as detailed in Sec. II A 4. The number of electronic structure calls for dTSH (and identically TSH) is constant at 435, as each of the 87 initial conditions was run five times. This number is comparable to the one obtained for ESSAIMS, even if it fluctuates to higher values when spawning events take place. In contrast, the average number of electronic structure calls for AIMS is significantly smaller, as there is only one run per initial condition. The AIMS dynamics starts with 87 necessary calls, and this number increases to 163 within 50 fs, before the simulation stops due to the instabilities discussed above. As a curiosity, we also report the theoretical cost of an AIMS calculation without the IFGA, that is, where all parents would be coupled from t = 0 (gray line, upper panel of Fig. 4). The advantage of applying the IFGA in multiple spawning simulations is striking: Even for such a trivial nonadiabatic process, the number of electronic structure calls without the IFGA would make an AIMS dynamics intractable. At t = 0, 3828 electronic structure calls per time step would be required, increasing to 7381 calls per time step after the first 50 fs.
The use of five runs per initial condition is meant to carefully converge the respective stochastic algorithm of the (E/O)SSAIMS of Chemical Physics ARTICLE scitation.org/journal/jcp or dTSH dynamics, but it naturally leads to an increase in computational effort-an effort above that of AIMS for the first part of the dynamics. Interestingly, using a single run per initial condition for dTSH ["dTSH * (1 run)" in Fig. 4] only leads to a minor alteration of the population trace. The ESSAIMS result obtained with only one run is similar for the beginning of the dynamics to the converged ESSAIMS one but deviates slightly after around 80 fs of simulation.

B. Fulvene-Nonadiabatic dynamics through a sloped conical intersection
Besides the challenges related to its electronic structure, the photodynamics of cyclopropanone discussed above is relatively straightforward and similarly captured by both (E/O)SSAIMS and (d)TSH. In contrast, fulvene exhibits a very stable electronic structure, but a rather complex nuclear dynamics. Upon excitation to its first electronic state S 1 , fulvene can encounter two different conical intersections: 31-35 a peaked conical intersection leading to an efficient deactivation to S 0 mediated by a twist of the C= =CH 2 moiety or a strongly sloped conical intersection reached by a stretching of the C= =CH 2 bond, leading to a transfer of the molecule to S 0 but also a possible reflection process bringing the molecule back to the S 1 state [see Fig. 1(b) for a schematic representation of the sloped intersection]. In the following, we will focus on the dynamics through the sloped intersection as the outcome of this process appears to be sensitive to the method employed, constituting a severe test for (E/O)SSAIMS (see Ref. 35 for more information).
The bottom panel of Fig. 5 shows the S 1 population traces for AIMS, ESSAIMS, and OSSAIMS, as well as TSH and dTSH. We note that additional (E/O)SSAIMS dynamics were conducted to select an optimal selection criterion (see the supplementary material). Despite the complexity of the reflection process, both ESSAIMS (ϵ = 10 −5 a.u.) and OSSAIMS (ϵ = 10 −2 ) accurately reproduce the initial decay of population and the first revival of S 1 population obtained with AIMS. OSSAIMS predicts the same amount of reflected population as AIMS, while ESSAIMS only slightly underestimates it. The stochastic-selection strategies do not fully capture the second, much weaker reflection process after 35 fs of dynamics using their respective selection criterion. Looking at the average number of TBFs during the dynamics, one can deduce that the stochastic selection algorithm only takes effect after 10 fs of dynamics, when 2 TBFs are present on average for all methods. Subsequently, the average number of TBFs in AIMS grows significantly, up to almost 5, while it decreases in ESSAIMS and OSSAIMS and remains well below 2. In contrast with the agreement between (E/O)SSAIMS and AIMS, the population trace predicted by (d)TSH differs significantly from that obtained with AIMS, with more than twice the population appearing in S 1 after the reflection process. [We note that the simulation parameters of dTSH can have an important influence on the population decay in the dynamics of fulvene (see Ref. 35

for more information).]
Matching the standard error of (E/O)SSAIMS with that of (d)TSH requires five runs for the former and seven for the latter for each initial condition. Comparing the theoretical number of electronic structure calls per time step for the different methods (middle panel of Fig. 5) reveals that (E/O)SSAIMS is computationally less expensive than (d)TSH, thanks to the lower number of runs required. Interestingly, AIMS is the least expensive method during the first half of the dynamics, until 25 fs of dynamics, at which points its computational effort rises above the one of the fully converged (E/O)SSAIMS. The number of electronic structure calls per time step in AIMS rises from 18 at the very beginning of the photodynamics to 298 after 42 fs.
The ESSAIMS population trace obtained with only one run per initial condition already agrees closely with the fully converged result-within the standard error of the fully converged ESSAIMS result for most of the simulation, except for the initial decay at 10 fs and during the short repopulation at 28 fs. Conversely, the dTSH dynamics with a single run shows significant deviations from its converged result, lying well outside the standard error for most of the dynamics. This example further highlights the difference between the stochastic processes in (E/O)SSAIMS and (d)TSH (as discussed in Sec. II A): In (d)TSH, the stochastic process is used to describe the nonadiabatic transitions per se, and its convergence is crucial in complex nonadiabatic processes like here with fulvene; in (E/O)SSAIMS, the stochastic processes mostly take place after the nuclear wavepacket branching following a conical intersection, while the nonadiabatic transition itself remains described at the AIMS level. of Chemical Physics ARTICLE scitation.org/journal/jcp

C. 1,2-dithiane-Numerous nonadiabatic transitions caused by nearly degenerate electronic states
The interesting nonadiabatic dynamics of 1,2-dithiane has been revealed in a study employing dTSH: 36 Upon photoexcitation, dithiane commences an ultrafast ring-opening process in S 1 mediated by the breaking of its S-S bond, which allows the molecule to reorganize and extend for some time until the S-S bond reforms within 300 fs (in line with earlier experimental work 74 ). This intriguing nuclear dynamics represents a challenge for nonadiabatic methods. Besides the evident challenge of describing ring-opening processes from an electronic structure perspective, the excited-state dynamics of dithiane leads to a situation where the three lowest electronic states can become nearly degenerate. The interplay between these electronic states and the nonadiabatic transitions resulting from their near-degeneracy requires a proper treatment of the coupled electron/nuclear dynamics. Hence, this molecule provides an ideal test for SSAIMS, as its nonadiabatic dynamics will produce many TBFs, and their inter-and intrastate interactions will be crucial for an accurate description of the electronic populations.
Following photoexcitation, AIMS predicts that the S 1 population begins to decay after 25 fs (bottom panel in Fig. 6), leading to an S 1 population drop to 40% within 60 fs. Subsequently, the S 1 population experiences a revival-up to 80% after 75 fs-and then decreases until around 40% where it stabilizes (with some oscillations). OSSAIMS with a high threshold of ϵ = 0.5 closely reproduces this behavior: The initial decay and reflection are adequately described with the only difference being that the S 1 population is overestimated during the revival process. After 100 fs, the S 1 population of OSSAIMS also oscillates around 40%, in close agreement to AIMS. The average number of TBFs (top panel in Fig. 6) rises almost exponentially for AIMS, reaching nearly 15 TBFs per initial condition after 155 fs of dynamics. For OSSAIMS with ϵ = 0.5, this number remains between 1 and 2 and does not surpass 2.5 TBFs. Furthermore, running only a single run per initial condition for OSSAIMS leads to an S 1 population decay in close agreement to our converged run employing five runs per initial condition (light red dashed curve in Fig. 6).
The decay of the S 1 population in dTSH starts earlier than in (OSS)AIMS, and the revival in the early part of the dynamics is not reproduced by the mixed quantum/classical method (green curve in Fig. 6). However, the S 1 population stabilizes at the same level as AIMS and OSSAIMS (40%) after 100 fs of dynamics. We note that the maximum standard error of the S 1 decay of OSSAIMS (using five runs per initial conditions) is reproduced by running each initial condition eight times with dTSH.
Interestingly, while AIMS, OSSAIMS, and dTSH depict a rather similar S 1 population decay, monitoring the S 0 population reveals larger deviations between the methods (Fig. 7). OSSAIMS and AIMS show that some population appears in S 0 after 50 fs of dynamics, rising to about 25% within 100 fs before plateauing between 20% and 25%. In contrast, dTSH predicts that the initial rise of the S 0 population takes place earlier (mirroring the behavior observed for the S 1 decay) and stabilizes at a higher population, ∼40% of population after 100 fs. ESSAIMS with a threshold of ϵ = 10 −3 a.u. appears to capture rather well the S 0 population dynamics up to 100 fs of dynamics (light blue dashed curve in Fig. 7). However, the S 0 population continues to grow until it plateaus at around 40% of population, similar to dTSH. The AIMS results could be recovered only by reducing the selection threshold of ESSAIMS to a value of ϵ = 10 −5 a.u. 75 OSSAIMS (ϵ = 0.5) drastically reduces the cost of the dynamics when a large number of TBFs are present after ∼200 fs of dynamics, while nevertheless reproducing the AIMS population dynamics correctly. The cost of ESSAIMS (ϵ = 10 −3 a.u.) is lower than that of OSSAIMS (ϵ = 0.5), and even below the theoretical cost of dTSH (as more runs are needed for the latter to achieve a similar level of convergence). However, decreasing the ESSAIMS threshold to ϵ = 10 −5 a.u. to reach AIMS accuracy leads to a dramatic increase in the computational cost, even higher than that of AIMS for a large part of the dynamics. Hence, OSSAIMS once again appears to offer a good compromise between accuracy and computational cost. It is interesting to note that the S 0 population trace obtained with ESSAIMS (ϵ = 10 −3 a.u.) starts to diverge from that of AIMS (at around 90 fs) shortly after the cost of AIMS greatly surpasses the cost of ESSAIMS (at around 75 fs). Adding to this observation the similarity between the dTSH and ESSAIMS (ϵ = 10 −3 a.u.) population trace, one can infer that the coupling between trajectories (absent in dTSH and limited in ESSAIMS with a high threshold) is likely to play a role in the last part of the dynamics. OSSAIMS, which adequately reproduces the AIMS S 0 population trace, includes more TBFs, and its theoretical cost is comparable to that of AIMS until the plateau is reached after 100 fs of dynamics.
The difference in performance between OSSAIMS and ESSAIMS observed above highlights the importance of the selection criterion. In OSSAIMS, the criterion is solely based on the overlap between TBFs. In ESSAIMS, the selection criterion depends on the off-diagonal elements of the Hamiltonian matrix (in the basis of TBFs). As such, the selection process in ESSAIMS depends on whether the TBFs under consideration evolve on the same state or different states-in the intrastate case, the Hamiltonian matrix element will contain the nuclear kinetic energy operator and the electronic energy, while in the interstate case the Hamiltonian matrix element contains the scalar product of the nonadiabatic coupling vectors with the nuclear gradient. In practice, we observe that ESSAIMS would be more likely to initiate a stochastic selection for TBFs evolving on different states than for TBFs on the same state, as the interstate coupling term is more likely to reach a small value (due to vanishing nonadiabatic coupling terms) than the intrastate one. OSSAIMS, on the other hand, does not differentiate between these two cases and only focuses on the overlap between TBFs. A closer look at the respective Hamiltonian matrix elements during the 1,2dithiane dynamics provides more insight into the difference between OSSAIMS and ESSAIMS. During the entire OSSAIMS dynamics, the overlap between TBFs remains rather large, and increasing the selection threshold to 0.7 does not significantly alter the dynamics (see the supplementary material). Hence, OSSAIMS will preserve any coupled TBFs, irrespective of their electronic state. Such couplings between TBFs appear to be critical to reproduce the AIMS dynamics. (We found that in order to recover the result of ESSAIMS with ϵ = 10 −3 a.u., it was necessary to increase the threshold of OSSAIMS to 0.8-see the supplementary material.) Focusing now on the offdiagonal elements of the Hamiltonian matrix, we observe that they range between 10 −3 a.u. and 10 −5 a.u. for most of the dynamics. The ESSAIMS thresholds used in our simulations thus represent two limiting cases-ϵ = 10 −3 a.u. is larger than the off-diagonal elements between any TBFs and will lead to an immediate stochastic selection, while ϵ = 10 −5 a.u. is a lower limit that the coupled TBFs only rarely reach. Interestingly, ESSAIMS cannot accurately reproduce the AIMS dynamics even with an intermediate threshold of ϵ = 10 −4 a.u. (see the supplementary material). This shows that the overlap between TBFs appears to be a robust criterion for SSAIMS, not only for the specific case of 1,2-dithiane but also for the other examples discussed above.

IV. SUMMARY AND CONCLUSION
In summary, we applied the novel framework of stochasticselection ab initio multiple spawning (SSAIMS) to the challenging photodynamics of different molecular systems to highlight its advantages and limitations and compare its performance with the celebrated mixed quantum/classical method TSH.
The results obtained for cyclopropanone indicate that both OSSAIMS and ESSAIMS can stabilize an AIMS dynamics suffering from electronic structure instabilities for weakly coupled TBFs. A very small selection threshold, i.e., a dynamics that remains very close to AIMS, could already achieve such a stabilization. By choosing an adequate selection threshold, the full decay of the S 1 population can be simulated with both flavors of SSAIMS. The computational cost of (E/O)SSAIMS remains close to that of TSH, with OSSAIMS requiring on average fewer TBFs for a result almost identical to ESSAIMS. The photodynamics of fulvene constitutes a stringent test for the robustness of (E/O)SSAIMS due to multiple passages through the S 1 /S 0 seam of intersection. Undeterred by such nonadiabatic processes, (E/O)SSAIMS could reproduce the AIMS population decay for the challenging photodynamics of fulvene, despite a rather aggressive selection of the TBFs. Moreover, (E/O)SSAIMS becomes rapidly cheaper than AIMS, even with numerous runs per initial conditions to ensure a full convergence of the stochastic algorithm. In contrast to (d)TSH, an even cheaper version of (E/O)SSAIMS, using a single run per initial conditions, still achieves an overall good agreement with the AIMS reference. The challenging photodynamics of 1,2-dithiane requires a large number of TBFs for its depiction. OSSAIMS was able to reproduce the dynamics predicted by AIMS while decreasing the computational cost significantly. Interestingly, ESSAIMS with a loose selection criterion deviates from AIMS and reproduces the dTSH dynamics at longer times, exhibiting the effect of mimicking the independent trajectory approximation of (d)TSH by considering the TBFs mostly as uncoupled.
Overall, both OSSAIMS and ESSAIMS proved to be stable and reliable strategies that, in many cases, could provide AIMS-quality nonadiabatic dynamics at a much-reduced cost, often competitive with the mixed quantum/classical method TSH. We showed that ESSAIMS and OSSAIMS achieve a very similar accuracy for both cyclopropanone and fulvene, with OSSAIMS necessitating slightly fewer TBFs on average while providing a slightly better agreement with AIMS. For 1,2-dithiane, OSSAIMS predicts the same dynamics as AIMS while reducing the necessary TBFs to below 2.5 on average. In comparison, ESSAIMS captures the early behavior of the dynamics well but, at later times, collapses to the dTSH result when using of Chemical Physics ARTICLE scitation.org/journal/jcp a threshold leading to a comparable cost with OSSAIMS. Achieving convergence to the AIMS result requires using a dramatically smaller selection criterion for ESSAIMS, resulting in an overall cost surpassing that of AIMS. Hence, these three exemplary molecular test systems indicate that OSSAIMS provides a more reliable and cost-efficient framework for further applications. It remains to be noted that both SSAIMS strategies require several test runs to determine an adequate stochastic-selection criterion. The development of an alternative, parameter-free version of SSAIMS is currently ongoing.

SUPPLEMENTARY MATERIAL
See the supplementary material for a validation of the active space chosen for 1,2-dithiane and additional (E/O)SSAIMS dynamics for fulvene and 1,2-dithiane with other selection criteria. All the initial conditions employed in this work as well as all the population traces are also provided as the supplementary material.