SSAIMSStochastic-Selection Ab Initio Multiple Spawning for Efficient Nonadiabatic Molecular Dynamics

Ab initio multiple spawning provides a powerful and accurate way of describing the excited-state dynamics of molecular systems, whose strength resides in the proper description of coherence effects during nonadiabatic processes thanks to the coupling of trajectory basis functions. However, the simultaneous propagation of a large number of trajectory basis functions can be numerically inconvenient. We propose here an elegant and simple solution to this issue, which consists of (i) detecting uncoupled groups of coupled trajectory basis functions and (ii) selecting stochastically one of these groups to continue the ab initio multiple spawning dynamics. We show that this procedure can reproduce the results of full ab initio multiple spawning dynamics in cases where the uncoupled groups of trajectory basis functions stay uncoupled throughout the dynamics (which is often the case in high-dimensional problems). We present and discuss the aforementioned idea in detail and provide simple numerical applications on indole, ethylene, and protonated formaldimine, highlighting the potential of stochastic-selection ab initio multiple spawning.


I. INTRODUCTION
Developing accurate and efficient nonadiabatic ab initio molecular dynamics methods is an important challenge for theoretical chemistry. At least two obstacles can be identified. First, the description of nonadiabatic phenomena requires nuclear dynamics beyond the Born−Oppenheimer approximation and, quite often, beyond the standard classical approximation for the nuclei. Second, the nuclear dynamics relies on accurate and efficient electronic structure calculations for potential energy surfaces, their gradients, and also derivatives of the electronic wave functions (which lead to the nonadiabatic coupling matrix elements that govern the propensity for transitions between electronic states). In this work, we focus primarily on the first problem, i.e. the efficient and accurate description of nonadiabatic dynamics. In principle, the required potential energy surfaces and couplings might come from parametrized analytic functions or "on the fly" ab initio quantum chemical methods, although we will generally prefer the latter.
Nonadiabatic dynamics methods are usually classified in different families. Grid-based wavepacket dynamics can provide a numerically exact description of nonadiabatic events, but unfortunately turns out to be computationally practical only for few nuclear degrees of freedom. 1,2 Trajectory-based methods, such as Ehrenfest dynamics or trajectory surface hopping 3−7 (TSH), propose a mixed quantum/classical description of nonadiabatic dynamics (treating the nuclear coordinates classically and providing a quantum mechanical description for the electronic coordinates or some coarse-grained representation thereof). Although the dynamics of molecules in their full dimensionality is rendered tractable by trajectory methods, this comes at the cost of uncontrolled (and generally difficult to improve) approximations.
At the intersection between these quantum grid and mixed quantum/classical trajectory methods lies the ab initio multiple spawning (AIMS) approach. In AIMS, the nuclear wave function is represented by a swarm of coupled frozen Gaussian functions following classical trajectories, whose evolution is computed "on-the-fly" with a given electronic structure method. 8−11 The number of trajectory basis functions (TBFs) used to describe the wavepacket increases during the dynamics (through "spawning") to accurately describe wavepacket bifurcation in nonadiabatic regions. In the limit of a complete basis set and exact evaluation of all the integrals required for the propagation, AIMS constitutes a formally exact solution of the time-dependent Schrodinger equation. Despite these attractive features, the number of coupled basis functions spawned during dynamics can become uncomfortably large. This is problematic because most of the computational burden in AIMS is due to electronic structure calculations, and the number of required electronic structure calculations grows at least linearly and at worst quadratically with the number of TBFs. 12 In many problems, nonadiabatic transitions are relatively rare events, and in these cases the increase in the number of TBFs causes little difficulty. However, there are important problems that involve a plethora of nonadiabatic transitions, especially when the density of electronic states is large. A prominent example is electronic energy transfer in multichromophoric systems such as photosynthetic light harvesting proteins 13−15 or conjugated polymers. 16 −18 In these problems, a myriad of electronic state crossings is expected and this will rapidly lead to a large number of TBFs in AIMS (see ref 19 for a recent strategy using TBFs for dense manifolds of electronic states). In some cases, population transfer may be highly efficient at these crossings, leading to the frustrating situation where there is an ever increasing number of TBFs, many of which carry little or no population (and thus add nothing to the physical description).
In this article, we propose a variation of the AIMS method, called stochastic-selection AIMS (SSAIMS), which exploits the decoupling between TBFs that occurs naturally in the course of the dynamics to reduce the number of running trajectories. SSAIMS still preserves the accurate description of nonadiabatic events provided by AIMS and can be converged to the AIMS results in well-defined limits. Numerical tests on different molecules are further presented to validate the SSAIMS method.
into the time-dependent Schrodinger equation leads to an equation of motion for a nuclear wave function Ω I (R,t) in In eq 2, Φ I (r; R) is the electronic wave function for state I with corresponding electronic energy E el I (R), TN is the nuclear kinetic energy operator, N is the number of nuclei, and ρ indexes the Cartesian coordinates of the nuclei. The electronic states are chosen to be orthonormal in eq 2 and the following.
Full multiple spawning (FMS) represents the nuclear wave function for each electronic state I by a linear combination of multidimensional frozen Gaussians called trajectory basis functions (TBFs): 11,21,22 A particular frozen Gaussian l evolving on electronic state I at time t, χ l I (R; R̅ l I (t), P̅ l I (t), γ ̅ l I (t), α), has its phase-space center located at position R̅ l I (t) and momentum P̅ l I (t). Its width α is time-independent (i.e., frozen) and the phase γ ̅ l I (t) is integrated semiclassically in order to factor out the largest oscillations from the complex amplitudes, C l I (t), and improve the numerical stability of their propagation. 10 A set of equations of motion for the complex amplitudes C l I (t) is obtained by inserting the FMS Ansatz into eq 2: The TBF centers simply follow classical trajectories on a particular electronic state, although other equations of motion are also possible (such as the Ehrenfest-type equations used in ab initio multiple cloning 23,24 ). In eq 4, bold symbols indicate matrices with respect to the nonorthogonal Gaussian basis functions. The overlap matrix S has elements (S II ) kl = ⟨χ k I |χ l I ⟩ R and its corresponding right-acting time-derivative Ṡis defined by In FMS, the Hamiltonian matrix mediate interstate and intrastate amplitude transfer, where the last two terms in eq 6 are responsible for nonadiabatic effects: We note that the second derivative G IJ terms are usually neglected, 10 as justified recently. 25 Usually, eq 4 is rewritten in terms of differences of electronic energies (with respect to a ground-state reference) in order to make the numerical integration more stable. A similar procedure is also often used in TSH. 26 A key feature of the FMS method is that it uses an adaptive basis set to ensure an accurate description of nonadiabatic processes. The number of TBFs describing the nuclear wave function for the Ith state, N I (t), changes in time as a result of spawning events. In short, a TBF entering a region of strong nonadiabaticity can spawn one or more new TBFs onto the coupled electronic state. Upon spawning, the size of the matrices in eq 4 is extended and the resulting coupled propagation of the enlarged set of TBFs allows for exchange of nuclear amplitude between electronic states. For detailed discussions about the spawning algorithm, the reader is referred to previous work. 8,10,27 The Journal of Physical Chemistry A pubs.acs.org/JPCA Article In AIMS, 8−10 eq 4 is solved along classical trajectories computed on-the-fly using a given electronic structure method, which provides the electronic energies and nonadiabatic coupling vectors required by eq 6. Two important numerical approximations simplify this coupled electron−nuclear propagation for molecular systems. First, the integrals in eq 6 are computed in a zeroth-order saddle-point approximation, 10 i.e., where R kl (c) is the position of the centroid between TBFs k and l (i.e., the position of the maximum of the product of the two TBFs). Second, the initial nuclear wave function is represented at time t = 0 by a given number of initial TBFs, which will be run independently in the subsequent nuclear dynamics (the independent first generation approximation). A number of reviews and previous papers provide extended discussions about both approximations. 10,28−33 II.b. Stochastic-Selection Ab Initio Multiple Spawning. At the beginning of an AIMS run, the initial TBF evolves on a given electronic state I. When this TBF encounters a region of strong nonadiabatic coupling with a different electronic state J, a new TBF is created on the Jth electronic state in a process called spawning. The two TBFs are now connected by an off-diagonal element of the H matrix in eq 6 (involving both nuclear and electronic coordinates, i.e. this is the vibronic matrix element H 12 IJ ). Each TBF can spawn additional TBFs during the dynamics, resulting in a growing number of coupled trajectories for a particular initial condition (see Figure 1). A large number of TBFs obviously implies a large number of electronic structure calculations that might complicate the application of AIMS to long-time dynamics, or nonadiabatic dynamics with a large number of electronic states. A further practical consideration is the possibility that one of the electronic structure calculations might fail for some numerical reason (e.g., divergence of iterative orbital optimization). This is a well-known difficulty in even ground-state ab initio molecular dynamics methods but tends to be more problematic when excited electronic states are involved, because of the added complexity when multiple electronic states need to be obtained at a given molecular geometry. This difficulty is further exacerbated when multiple TBFs are being propagated, because a failure in any of the electronic structure calculations will necessitate stopping the entire dynamics simulation and intervention (ideally automated, but often manual) to solve the problem. The probability of requiring such intervention will therefore increase with the number of TBFs being propagated.
In order to reduce the number of costly electronic structure calculations and simultaneously minimize the probability of manual intervention, it would be highly desirable to reduce the number of TBFs during the simulation. This could be accomplished with a "death" process that decreases the basis set size, in analogy to the "spawning" process that increases the basis set size. Indeed, such death processes have been used previously, but normally only when the dynamics is deemed uninteresting. For example, when one is only interested in the excited-state lifetime, there is rarely any need to follow groundstate TBFs for the entire length of the simulation. 34,35 Here, we introduce a general stochastic death process that in many cases does not change observable results at all. The key observation is that the expectation value of an operator Ôfor the nuclear wave function associated with the Ith electronic state is calculated in AIMS as Figure 1. Schematic representation of a typical SSAIMS propagation for a particular initial TBF (red Gaussian function with center depicted by the red circle). Initially, a single TBF evolves on a given electronic state (thin red line). A spawn occurs after a short time and two TBFs are now propagated, coupled by the Hamiltonian matrix given by eq 6. A horizontal dashed line represents the coupling between two TBFs (the stronger the coupling, the darker the line). Additional TBFs are spawned over the course of the dynamics. Until t SS , the SSAIMS algorithm does not alter the dynamics in any way and all TBFs are directly or indirectly coupled. A detailed pictorial representation of the stochastic-selection process at time t SS is presented on the right. (a) The three TBFs on the left are coupled, while the one on the right is only weakly coupled to the third TBF. (b) The stochastic-selection algorithm detects the weak coupling (coupling is below the threshold ε) and forms two uncoupled blocks of coupled trajectories (3 TBFs in the first block, 1 in the other). (c) A random number is generated and one of the blocks is selected (with selection probability proportional to the population in the block). In this example, the first block is selected. The unselected blocks (the fourth TBF in this example) are then removed from the simulation and the population of each TBF in the selected block is scaled to maintain normalization of the wave function. The nonadiabatic dynamics then proceeds with only TBFs from the selected block.
Most observables of interest correspond to local operators, implying that the matrix elements in the numerator will vanish when the bra and ket TBFs are well-separated in phase space. Thus, one only needs to consider pairs of TBFs that are close to each other. Typically, this is accomplished by "screening," which neglects the observable matrix elements in the numerator when the corresponding overlap matrix element is small. In the limit where all TBFs are distant from each other, eq 11 reduces to a fully incoherent ("classical-like") sum: where we have defined the population of a basis function as Similar considerations apply to the nonadiabatic coupling matrix elements, i.e. two TBFs will be most effectively coupled when they are close to each other (implying that the overlap matrix element for the pair of TBFs is non-negligible) and also close to an avoided crossing or conical intersection (where the nonadiabatic coupling is large). These observations suggest a new approach to AIMS dynamics that incorporates sparsity: whenever TBFs (or groups of TBFs) become uncoupled, the full AIMS dynamics becomes superfluous as there is no population transfer or observable quantum interference between the uncoupled groupsthis observation is reminiscent of coarse graining in quantum decoherence introduced by Rossky and co-workers. 36 Analysis of eq 11 suggests that one could split the "mother" simulation into several "daughter" simulations (one for each group of uncoupled TBFs). As long as the uncoupled groups do not become coupled later in the dynamics, any observables computed using eq 11 for the set of daughter simulations (summing the results across simulations with appropriate weights) will be identical to the results that would have been obtained without splitting. This is shown schematically in Figure  2, where the dotted lines indicate that the dynamics from the "mother" simulation is being reused without recomputation. The daughter simulations could be run in parallel given sufficient computational resources, or they could be run sequentially. The total number of TBFs in the simulation will not change compared to a full AIMS run. However, this splitting approach minimizes the number of TBFs in any individual simulation, which ensures that (1) a failure in any one of the TBF electronic structure calculations affects as few other TBFs as possible, (2) the dimension of the matrices is minimized (simplifying inversion of the overlap matrix), and (3) the number of required off-diagonal Hamiltonian matrix elements is minimized. In practice, the benefits of the latter two considerations are less than might have been expected from a formal analysis, since AIMS already uses sparsity (through a conjugate gradient procedure) to solve eq 4 and also avoids computing negligible Hamiltonian matrix elements by the screening procedure described above.
We will christen the splitting procedure just described as "AIMS with splitting." Simulations are split into uncoupled groups of TBFs (as these are identified), and then the "daughter" simulations are run independently. Each of the daughter simulations is split into uncoupled groups of TBFs, leading to "granddaughter" simulations that are run independently. The process repeats and all simulations are collected together at the endtheir union is an approximation to the full AIMS dynamics that would have been obtained without splitting. Although this can already be advantageous, the procedure does have the drawback of considerable bookkeeping since the initial conditions for each of the daughter simulations need to be recorded and ultimately one needs to go back and finish the propagation for each of these (or use increasing computational resources throughout the simulation to run these in parallel). A stochastic approach may be more convenient. In this case, one randomly chooses (with appropriate probability) one of the daughter simulations at each splitting event until the end of the propagation time. Repeating this procedure with the same initial conditions many times and taking the union of the resulting TBFs (appropriately weighted) would lead to the same result as obtained in AIMS with splitting. The drawback here is that numerous paths may be replicated, but the advantage is that there is no bookkeeping. We call this "SSAIMS with repetition." Figure 3 schematically compares these different approaches. Both AIMS with splitting and SSAIMS with repetition ensure that the fully coupled AIMS propagation is used only when needed, thus reducing the number of evolving TBFs in a given run. Finally, we note that one can also independently stochastically sample both initial conditions and the choice of uncoupled groups at each splitting event. Although this is the variant that we will recommend (SSAIMS with repetition, independently sampling initial conditions and selection of TBF groups), we present results for some of the other variants to fully characterize convergence.
The idea of stochastic selection rests on the approximation that by stochastically selecting one TBF, or a group of TBFs, we consider that it will remain uncoupled to the rest of the AIMS swarm for all future times. In other words, we consider that uncoupled TBFs will never overlap again during the nonadiabatic dynamics. This "no-reoverlap" approximation is closely related to the notion of Poincarérecurrence time. We Figure 2. Schematic splitting procedure for AIMS simulations. The "mother" simulation proceeds until the TBFs are naturally arranged in weakly coupled groups. At this point, each of the TBF groups forms the initial conditions for an independent "daughter" AIMS dynamics simulation. For each of the daughter simulations, the TBFs again separate and the simulation can be split into "granddaughter" simulations. Dotted lines indicate parts of the daughter simulation that are not replicated as they are present in the mother simulation. In many multidimensional problems, TBFs will not come close to each other after separating. In that case, the union of all the TBFs created in the mother/daughter/granddaughter simulations will replicate the results that would have been obtained without simulation splitting.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article expect it to be an excellent approximation for high-dimensional systems, but is likely to fail for low-dimensional ones. Let us now define more precisely what we mean by uncoupled. In the following, we will compare the absolute value of the matrix elements H kl ·· (both H kl IJ and H kl II ) with a threshold criterion to determine whether trajectories can be considered uncoupled. Hence, the stochastic selection will be applied on truly uncoupled TBFs only if a sufficiently small threshold is selected. On the other hand, we can relax the definition of uncoupled TBFs and perform the stochastic selection more often by employing a larger threshold. The resulting reduced number of TBFs in the dynamics comes at the cost of neglecting the effect of nonvanishing interference terms. Since the proposed criterion has an energy unit, we label the resulting stochastic-selection AIMS method as ESSAIMS. We finally note that other definitions of uncoupled TBFs could obviously be used for SSAIMS, such as the overlap between two TBFs. This definition was tested (see the Supporting Information) and found to closely reproduce ESSAIMS results.
The ESSAIMS method works as follows. After each nuclear time step, the algorithm monitors the off-diagonal matrix elements of H between the running TBFs (heavy dashed lines in Figure 1). If the absolute value of one matrix element falls below a preset threshold ε, the corresponding TBFs are considered uncoupled (Figure 1a). 37 More precisely, a connectivity matrix M between TBFs is constructed and element (M) kl is set to unity if |(H) kl | > ε, i.e., if TBFs k and l are coupled according to the predefined threshold, and zero otherwise. The algorithm should also detect indirect couplings, which are couplings between two TBFs occurring through one or more other TBFs. These couplings are easily identified by constructing powers of the connectivity matrix. If TBFs k and l are indirectly coupled via one other TBF m, the matrix element (M 2 ) kl will be nonzero. If the indirect coupling between k and l takes place through two other TBFs, the matrix element (M 3 ) kl will be nonzero. These relationships continue for indirect couplings of higher orders. The matrix multiplication is therefore performed until no more changes happen in the connectivity pattern. Special care is needed during a spawning process, as coupled propagation with the newly created TBF is usually initiated in a region where couplings are small. The stochastic-selection algorithm is therefore suspended from the entry time of a new TBF until spawning time, preventing the newly created TBF from generating an unwanted selection event. It is finally important to remark that the ε threshold is likely to be system dependent (since it is given in energy units), and will need to be determined prior to any ESSAIMS dynamics. It is however clear that the asymptote for the off-diagonal Hamiltonian matrix elements will always tend toward zero whenever TBFs are sufficiently distant.
Once the coupling pattern between TBFs is identified, uncoupled blocks of coupled trajectories are identified (Figure 1b) where ζ is a random number generated in the interval [0: 1]. All TBFs that are not in the βth block are removed from the simulation and the coefficients of the TBFs in the βth block are scaled to maintain wave function normalization, as depicted in Figure 1c. The simulation then continues with the reduced subspace of TBFs. The ESSAIMS method therefore ensures a reduced number of TBFs throughout the dynamics for any given simulation when a sufficiently large value of ε is used. ESSAIMS can either be used as an approximation to AIMS by running only a few stochastic trajectories, or it can be fully converged to recover the AIMS population distribution among the different electronic states. The first case is particularly attractive for the nonadiabatic dynamics of large molecular systems or for dynamics involving several excited states. The second case would be useful if important electronic structure instabilities are encountered during an AIMS dynamics. If the different uncoupled blocks of TBFs are run independently to completion (AIMS with splitting), the scheme is reminiscent of the "ant" method in surface hopping. 3, 38 We note here that another approach, coined apoptosis, was recently proposed to remove moving nonorthogonal basis functions and stabilize matrix inversion in quantum dynamics simulations. 39 It is finally interesting to note that ESSAIMS shares some similarities with the TSH algorithm commonly used in nonadiabatic dynamics. In TSH, the nuclear dynamics is represented by a swarm of uncoupled classical trajectories that can hop from one electronic state to another in nonadiabatic regions. The hopping process is determined by means of a stochastic algorithm, based on probabilities computed after each nuclear time step from propagated amplitudes in different electronic states. As the trajectories are all independent, any expectation value is obtained by definition from an incoherent analysis. In some sense, the selection of ESSAIMS resembles the hopping process of TSH, as a swarm of TBFs evolving in different states are suddenly reduced to a subset of renormalized TBFs whenever the subset can be considered independent from the other TBFs. Nevertheless, two major differences should be highlighted. First, the ESSAIMS calculation maintains coupling between TBFs as they separate−this captures decoherence effects that are often approximated (and sometimes ignored) in TSH. Second, a fully converged ESSAIMS calculation (small The Journal of Physical Chemistry A pubs.acs.org/JPCA Article coupling threshold ε) will reproduce the AIMS results. With sufficient fully coupled initial conditions, this is guaranteed to converge to exact quantum dynamics (although it may be difficult to reach this limit in practice). In TSH, decoherence effects are only captured with ex post facto corrections. 36,43−48 The ambiguity in describing decoherence in TSH stems from the lack of a rigorous derivation for the procedure, although there have been some notable attempts. 40−42 Furthermore, TSH is not guaranteed to converge to the exact quantum dynamical result, regardless of the number of trajectories (one consequence of this is the well-known dependence of TSH results on the electronic representation used, i.e. adiabatic or diabatic). II.c. Computational Details. The nonadiabatic dynamics of ethylene, protonated formaldimine, and indole were used to test the ESSAIMS method and to compare it with AIMS. AIMS and ESSAIMS dynamics for the ethylene and protonated formaldimine were performed with a modified version of the AIMS/ MolPro 2012.1 interface. 12,49 The time step used for the classical propagation was 20 atomic time units (atu) (reduced to 5 atu in coupling regions). The spawning threshold was fixed to 3.0 a.u. −1 (magnitude of the nonadiabatic coupling vectors) and the minimum population to spawn was set to 0.01 (0.05 for protonated formaldimine). Initial conditions (positions and momenta representing the center of the initial TBFs) were randomly sampled from a harmonic Wigner distribution, with 60 samples for ethylene and 39 for protonated formaldimine. The electronic structure of ethylene was described with SA-3-CASSCF(2/2)/6-31G*, i.e. state average complete active space self-consistent field (SA-CASSCF) calculation with an active space of two electrons in two orbitals, equally weighting the lowest three singlet states (S 0 , S 1 , and S 2 ) and using the 6-31G* basis set. 50−52 For the protonated formaldimine molecule, SA-4-CASSCF(6/5)/6-31G* was used.
The AIMS dynamics of indole was computed with the AIMS/ TeraChem interface. 53 Classical propagation used a time step of 20 atu (reduced to 5 atu in coupling regions). The spawning threshold was set to 0.005 au (scalar product of the nonadiabatic coupling vector and velocity), and the minimum required population to spawn was 0.1. The first three electronic states were considered in the dynamics and the electronic structure was calculated with density functional theory 54,55 (DFT) and linear-response time-dependent density functional theory 56−58 (LR-TDDFT), applying the Tamm−Dancoff 59,60 and adiabatic approximations. We used the PBE 61 exchange and correlation functional and the 6-311G 62 basis set. As a consequence of the adiabatic approximation, 63 LR-TDDFT poorly describes the coupling between the ground state and the first excited electronic state. 64,65 This problem, however, is not relevant here since we focus on transitions from S 2 to S 1 (and not subsequent relaxation to S 0 ). The initial conditions for indole were obtained by first equilibrating the system for 7.9 ps with NVT (300 K) ground-state ab initio molecular dynamics in TeraChem 66−68 using DFT-PBE/6-311G. A frame was then extracted from the end of a subsequent NVE run (18.4 ps) at the same level of theory. All error bars indicate standard errors of the mean.

III. RESULTS
We validate the ESSAIMS approach by comparing to AIMS for excited-state dynamics of various molecules. We start by discussing the difference between an AIMS and an ESSAIMS run for the excited-state dynamics of a medium-size molecule, indole. We then address the question of ESSAIMS performance based on the ethylene and protonated formaldimine nonadiabatic dynamics. First, we determine how many ESSAIMS runs are necessary to reproduce the AIMS population trace. Each of the ESSAIMS runs is initiated from the same initial condition (positions and momenta), but with a different seed for the random number generator (which influences the stochasticselection process). Second, we study the influence of the ε threshold parameter on the population trace for a fixed number of runs. Finally, we compare the average ESSAIMS population trace with AIMS, by using a small number of ESSAIMS runs for each initial condition.
It is worth noting at this stage that the purpose of this study is to validate the ESSAIMS methodology and not to provide a full We start by illustrating the behavior of SSAIMS with the excited-state dynamics of indole. In the following, we focus on the dynamics produced by SSAIMS for a single initial condition and compare it with AIMS (as done schematically in Figure 3). This first example is meant to illustrate the stochastic-selection process and exemplify some potential limitations of the method. It is important to note at this stage that a single SSAIMS run is not expected to reproduce the AIMS dynamics for a single initial condition, and multiple runs are required as discussed above. Numerical examples showing how SSAIMS can be converged toward AIMS will be presented later.
The AIMS dynamics of indole along a given initial condition is presented in Figure 4 (left panel). The population of the S 2 state decreases rapidly after less than 500 atomic time units (atu) of dynamics but reappears briefly before a final transfer to S 1 occurs at the end of the run. A specific ESSAIMS run also predicts a complete depopulation of the S 2 state by the end of the run (Figure 4 left panel), but the mechanism varies depending on the threshold used (and will also vary depending on the random seed). For a better visualization of the selection events, we use a projection of the running TBFs on the central CC bond length of indole (right panel, Figure 4). Each line represents a TBF and the thickness of the line is directly proportional to the population of the TBF. For a medium threshold (10 −4 a.u.), the stochastic-selection algorithm enters into play directly after the first nonadiabatic event (before 500 atu). In this example, the algorithm dismisses the TBF running on S 2 and only a single TBF evolves on S 1 until the end of the run, with a population renormalized to 1.0. For a slightly larger threshold parameter (10 −5 a.u., middle panel on the right of Figure 2), the stochastic selection algorithm does not consider the two TBFs as uncoupled after the first transfer, and a second exchange of population (without any new spawning event) takes place at 700 atu. As a result, the second TBF transfers all its population back to the first TBF on S 1 . After 1250 atu of dynamics, a new spawn occurs from the first TBF and a complete transfer to S 1 takes place. The selection algorithm then detects that all TBFs in S 1 are coupled together, but none is coupled with the one in S 2 , resulting in a selection of the TBFs in S 1 (the line width slightly changes in Figure 4 as a result of the renormalization). Although the overall population transfer is well captured by ESSAIMS with a 10 −5 threshold, it misses the last spawning event taking place 1800 atu in the AIMS dynamics (upper panel on the right of Figure 4), leading to a complete depopulation of S 2 . Importantly, changing the seed for the random number generator can alter the overall dynamics, as a different group of coupled TBFs might be selected in the stochastic process. Therefore, convergence of ESSAIMS toward the AIMS results requires an average over multiple runs, each one being initiated with a different seed (see following sections).
In terms of performance, the AIMS dynamics required 1619 electronic structure evaluations (12.2 h on 2 GeForce GTX Titan GPU cards), while the ESSAIMS dynamics only required 1211 electronic structure evaluations (9.0 h on the same machine) with ε = 10 −5 a.u. or 354 evaluations (2.5 h wall time) with ε = 10 −4 a.u.
Despite the fact that the overall depopulation of S 2 is well predicted by the ESSAIMS dynamics, the fine details of the nonadiabatic process will be lost when a single ESSAIMS run is compared to AIMS starting with the same initial conditions. The threshold value is also likely to be system dependent and initial tests are needed to estimate its value. A small threshold will obviously capture all interferences between TBFs (with a vanishing threshold reproducing AIMS results exactly) but will then require many TBFs to be treated simultaneously. On the other hand, a larger threshold will significantly decrease the computational cost but may considerably alter the nonadiabatic processes, as illustrated by the indole dynamics with ε = 10 −4 a.u. In this particular case, the stochastic selection is triggered too rapidly and does not allow for a second interaction between the two TBFs at a later time. The large threshold therefore compromises the accuracy of the no-reoverlap approximation discussed above. This example shows that one must carry out numerous samples when using ESSAIMS, and this is pursued in the following.
III.b. Nonadiabatic Decay of Ethylene and Protonated Formaldimine. The S 1 population decay of ethylene (H 2 C CH 2 ) is presented in Figure 5 (red curve) for a single initial condition in AIMS. We are first interested in reproducing this population trace by averaging over a number of ESSAIMS runs, initiated from the same initial condition but with different seeds for the random number generator. For a small enough ε threshold (10 −10 a.u.), only a few ESSAIMS runs are necessary to capture the overall dynamics of the S 1 population ( Figure 5, left panel). As expected, averaging over more runs converges to the AIMS result. The effect of the ε threshold for a fixed number of runs (50) is given in the right panel of Figure 5. Interestingly, the The Journal of Physical Chemistry A pubs.acs.org/JPCA Article overall averaged population varies only weakly for small values of ε, and sizable deviations (larger than a few percent) are observed only when the threshold is set to 10 −2 a.u. Based on these results, we computed an average of five runs for each initial condition used to obtain the total AIMS population trace ( Figure 6, red curves). The overall agreement between the fully coupled AIMS dynamics and ESSAIMS for the population decay of S 1 is excellent ( Figure 6, left panel), even for an aggressive threshold (10 −2 a.u.). During the nonadiabatic dynamics, the second excited state is slightly populated and ESSAIMS predicts this weak population transfer rather accurately ( Figure 6, right panel). It is finally interesting to compare the average number of TBFs per initial condition for each method (inset of the left panel, Figure 6). While AIMS requires an average of 8 TBFs (per initial condition) by the end of the dynamics, an average of fewer than 2 simultaneous TBFs are required using ESSAIMS. Furthermore, the average number of TBFs grows linearly with time in this AIMS example, while the average number of simultaneous TBFs required in ESSAIMS is roughly constant with propagation time. Since the number of electronic structure calculations required by AIMS is formally quadratic in the number of TBFs, ESSAIMS constitutes an appealing solution to reduce the cost of AIMS dynamics.
The accuracy of ESSAIMS in reproducing the AIMS result is further confirmed by the nonadiabatic dynamics of the protonated formaldimine (H 2 CNH 2 + ), initiated from its second electronic excited state. Upon photoexcitation, protonated formaldimine rapidly relaxes from S 2 to S 1 (red curve in Figure 7), losing 50% of S 2 population in the first 415 atu (10 fs) of dynamics, in agreement with former TSH results. 69,70 The ground state already starts to become populated after 460 atu (11 fs) of dynamics, but we will focus here on the S 2 and S 1 population traces (see the SI for detailed information about the populations of all electronic states).
As observed for the case of ethylene, ESSAIMS accurately mirrors the AIMS population trace for a single initial condition, even with rather high threshold values (see the SI). Furthermore, averaging over a small number of ESSAIMS runs (5 in average) with different seeds for each initial condition reproduces the S 2 and S 1 population trace given by AIMS, well within the error bars (Figure 7). In fact, a single run per initial condition is already sufficient to reproduce the overall transfer of population between S 1 and S 2 (blue curve in Figure 7). By limiting the average number of TBFs (inset in left panel of Figure 7), ESSAIMS reduces by a factor of 3 the average number of electronic structure calls per run with respect to AIMS: 1068 for AIMS (average over 39 initial conditions) compared to 316 for ESSAIMS (average over 194 runs, ε = 10 −3 a.u.).
In summary, the numerical tests presented here show that with appropriate numerical parameters (threshold and number of runs per initial condition), ESSAIMS accurately reproduces the AIMS population traces and allows for a substantial decrease in the number of coupled TBFs during the dynamics.

IV. CONCLUSIONS
The main conclusions of this work may be summarized as follows. We showed that a stochastic-selection algorithm can reduce the computational effort of AIMS dynamics by preventing the proliferation of weakly coupled TBFs during a simulation, leading to a stabilization of the average number of TBFs during a multiple spawning run. When nuclear basis functions can be grouped into blocks that are mutually uncoupled, the proposed algorithm stochastically picks one of the blocks to continue the nonadiabatic dynamics. In the limit of  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article truly uncoupled TBFs and long recurrence time, ESSAIMS is formally exact and can be converged to the AIMS result by performing several runs for each initial starting condition. We confirmed this by different numerical tests and showed that ESSAIMS reproduces accurately the AIMS results with a limited number of additional runs, while avoiding an exponential increase in the number of trajectories treated simultaneously, leading to a substantial reduction in the number of running TBFs. More generally, stochastic-selection AIMS opens new perspectives for the spawning process in nonadiabatic dynamics. While the spawning algorithm in AIMS tries to minimize the number of spawning events by judiciously placing the newly created child TBF, stochastic-selection AIMS allows for more general spawning criteria, as only importantor coupled TBFs will survive the selection process. Therefore, spawning could potentially take place as soon as a new TBF is expected to be useful, for not only nonadiabatic but also tunneling events, 73 since the selection process will discard it at a later time if this TBF is not important. This perspective will stimulate more work on stochastic-selection AIMS.
Additional information on the excited-state dynamics of the protonated formaldimine model, including the population trace obtained from AIMS, and an additional comparison between AIMS, ESSAIMS, and OSSAIMS (PDF)