Open Submitted: 27 January 2021 Accepted: 18 February 2021 Published Online: 09 March 2021
J. Chem. Phys. 154, 104110 (2021); https://doi.org/10.1063/5.0045572
more...View Affiliations
View Contributors
• Lea M. Ibele
• Yorick Lassmann
• Todd J. Martínez
• Basile F. E. Curchod

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 simultaneously 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 obtaining 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.
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,21. J. C. Tully, Theor. Chem. Acc. 103, 173 (2000). https://doi.org/10.1007/978-3-662-10421-7_32. M. Born and R. Oppenheimer, Ann. Phys. 389, 457 (1927). https://doi.org/10.1002/andp.19273892002 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 electron-nuclear 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–63. H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990). https://doi.org/10.1016/0009-2614(90)87014-i4. M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324, 1 (2000). https://doi.org/10.1016/S0370-1573(99)00047-25. G. A. Worth, H. D. Meyer, and L. S. Cederbaum, in Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, Advanced Series in Physical Chemistry Vol. 15, edited by W. Domcke, D. R. Yarkony, and H. Köppel (World Scientific, 2004), Chap. 14, pp. 583–617.6. H.-D. Meyer, F. Gatti, and G. A. Worth, Multidimensional Quantum Dynamics (John Wiley & Sons, 2009). 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 gold-standard 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. 77. G. W. Richings and S. Habershon, J. Chem. Phys. 148, 134116 (2018). https://doi.org/10.1063/1.5024869 and 88. G. W. Richings and S. Habershon, J. Phys. Chem. A 124, 9299 (2020). https://doi.org/10.1021/acs.jpca.0c06125 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,109. M. Barbatti, Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 620 (2011). https://doi.org/10.1002/wcms.6410. J. C. Tully, J. Chem. Phys. 137, 22A301 (2012). https://doi.org/10.1063/1.4757762 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–139. M. Barbatti, Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 620 (2011). https://doi.org/10.1002/wcms.6411. J. C. Tully, J. Chem. Phys. 93, 1061 (1990). https://doi.org/10.1063/1.45917012. S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys. 101, 4657 (1994). https://doi.org/10.1063/1.46745513. R. Crespo-Otero and M. Barbatti, Chem. Rev. 118, 7026 (2018). https://doi.org/10.1021/acs.chemrev.7b00577 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,1514. F. Agostini, S. K. Min, A. Abedi, and E. K. U. Gross, J. Chem. Theory Comput. 12, 2127 (2016). https://doi.org/10.1021/acs.jctc.5b0118015. S. K. Min, F. Agostini, I. Tavernelli, and E. K. U. Gross, J. Phys. Chem. Lett. 8, 3048 (2017). https://doi.org/10.1021/acs.jpclett.7b01249
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–1816. G. A. Worth, M. A. Robb, and I. Burghardt, Faraday Discuss. 127, 307 (2004). https://doi.org/10.1039/b314253a17. B. Lasorne, M. J. Bearpark, M. A. Robb, and G. A. Worth, Chem. Phys. Lett. 432, 604 (2006). https://doi.org/10.1016/j.cplett.2006.10.09918. G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, Int. Rev. Phys. Chem. 34, 269 (2015). https://doi.org/10.1080/0144235x.2015.1051354 multiconfiguration Ehrenfest (MCE),19–2119. D. V. Shalashilin, J. Chem. Phys. 130, 244101 (2009). https://doi.org/10.1063/1.315330220. K. Saita and D. V. Shalashilin, J. Chem. Phys. 137, 22A506 (2012). https://doi.org/10.1063/1.473431321. D. V. Makhov, C. Symonds, S. Fernandez-Alberti, and D. V. Shalashilin, Chem. Phys. 493, 200 (2017). https://doi.org/10.1016/j.chemphys.2017.04.003 and full and ab initio multiple spawning (FMS and AIMS).22–2422. T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100, 7884 (1996). https://doi.org/10.1021/jp953105a23. T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. A 101, 6389 (1997). https://doi.org/10.1021/jp970842t24. T. J. Martínez and R. D. Levine, J. Chem. Soc., Faraday Trans. 93, 941 (1997). https://doi.org/10.1039/a605958i 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.2525. B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martínez, Chem. Phys. 347, 3 (2008). https://doi.org/10.1016/j.chemphys.2008.01.014 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),2626. B. F. E. Curchod, W. J. Glover, and T. J. Martínez, J. Phys. Chem. A 124, 6133 (2020). https://doi.org/10.1021/acs.jpca.0c04113 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 excited-state 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,2-dithiane 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–3027. G. Cui, Y. Ai, and W. Fang, J. Phys. Chem. A 114, 730 (2010). https://doi.org/10.1021/jp908936u28. G. Cui and W. Fang, J. Phys. Chem. A 115, 1547 (2011). https://doi.org/10.1021/jp110632g29. M. Filatov, S. K. Min, and C. H. Choi, Phys. Chem. Chem. Phys. 21, 2489 (2019). https://doi.org/10.1039/c8cp07104g30. J. Suchan, J. Janoš, and P. Slavíček, J. Chem. Theory Comput. 16, 5809 (2020). https://doi.org/10.1021/acs.jctc.0c00512 Photoexcitation to the first excited singlet state (S1) triggers a ring-opening reaction, mediated by carbon–carbon bond breaking followed by subsequent dissociation into ethylene and CO. The S0 and S1 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 S1 to S0 without a significant contribution of population back transfer to S1. 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 S1 and S031–3431. M. J. Bearpark, F. Bernardi, M. Olivucci, M. A. Robb, and B. R. Smith, J. Am. Chem. Soc. 118, 5254 (1996). https://doi.org/10.1021/ja954279932. S. Alfalah, S. Belz, O. Deeb, M. Leibscher, J. Manz, and S. Zilberg, J. Chem. Phys. 130, 124318 (2009). https://doi.org/10.1063/1.308954633. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. J. Bearpark, and M. A. Robb, Phys. Chem. Chem. Phys. 12, 15725 (2010). https://doi.org/10.1039/c0cp01757d34. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. A. Robb, and M. J. Bearpark, J. Chem. Phys. 137, 22A548 (2012). https://doi.org/10.1063/1.4765087 and has subsequently been compared to a molecular version of Tully’s third model.3535. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f Figure 1(b) shows a one-dimensional cut of the S0 and S1 potential energy surfaces in the vicinity of the sloped conical intersection. Coming from the Franck–Condon region in S1, the nuclear wavepacket crosses the sloped conical intersection and, soon after being reflected on S0, splits into a part decaying to the S0 Franck–Condon region and a part transferring back to the S1 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 employed3535. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f 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.3636. C. D. Rankine, J. P. F. Nunes, M. S. Robinson, P. D. Lane, and D. A. Wann, Phys. Chem. Chem. Phys. 18, 27170 (2016). https://doi.org/10.1039/c6cp05518d By looking at the three lowest singlet states along a linear interpolation in internal coordinates between the Franck–Condon region, the S1 minimum geometry, and the S1/S0 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 ring-opening 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,3837. B. F. E. Curchod and T. J. Martínez, Chem. Rev. 118, 3305 (2018). https://doi.org/10.1021/acs.chemrev.7b0042338. L. M. Ibele, A. Nicolson, and B. F. E. Curchod, Mol. Phys. 118, e1665199 (2020). https://doi.org/10.1080/00268976.2019.1665199 TSH,13,3913. R. Crespo-Otero and M. Barbatti, Chem. Rev. 118, 7026 (2018). https://doi.org/10.1021/acs.chemrev.7b0057739. J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016). https://doi.org/10.1146/annurev-physchem-040215-112245 and nonadiabatic dynamics in general4040. F. Agostini and B. F. E. Curchod, Wiley Interdiscip. Rev. Comput. Mol. Sci. 9, e1417 (2019). https://doi.org/10.1002/wcms.1417 for additional details.
1. 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.1111. J. C. Tully, J. Chem. Phys. 93, 1061 (1990). https://doi.org/10.1063/1.459170 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–4841. E. R. Bittner and P. J. Rossky, J. Chem. Phys. 103, 8130 (1995). https://doi.org/10.1063/1.47017742. M. Thachuk, M. Y. Ivanov, and D. M. Wardlaw, J. Chem. Phys. 109, 5747 (1998). https://doi.org/10.1063/1.47719743. J.-Y. Fang and S. Hammes-Schiffer, J. Phys. Chem. A 103, 9399 (1999). https://doi.org/10.1021/jp991602b44. A. W. Jasper, S. Nangia, C. Zhu, and D. G. Truhlar, Acc. Chem. Res. 39, 101 (2006). https://doi.org/10.1021/ar040206v45. J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 024105 (2011). https://doi.org/10.1063/1.350677946. J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 244114 (2011). https://doi.org/10.1063/1.360344847. M. Persico and G. Granucci, Theor. Chem. Acc. 133, 1526 (2014). https://doi.org/10.1007/s00214-014-1526-148. B. F. E. Curchod and I. Tavernelli, J. Chem. Phys. 138, 184112 (2013). https://doi.org/10.1063/1.4803835 and a variety of ad hoc corrections have been proposed to alleviate this issue.39,45,46,49–5239. J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016). https://doi.org/10.1146/annurev-physchem-040215-11224545. J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 024105 (2011). https://doi.org/10.1063/1.350677946. J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 244114 (2011). https://doi.org/10.1063/1.360344849. G. Granucci, M. Persico, and A. Zoccante, J. Chem. Phys. 133, 134111 (2010). https://doi.org/10.1063/1.348900450. H. M. Jaeger, S. Fischer, and O. V. Prezhdo, J. Chem. Phys. 137, 22A545 (2012). https://doi.org/10.1063/1.475710051. N. Shenvi and W. Yang, J. Chem. Phys. 137, 22A528 (2012). https://doi.org/10.1063/1.474640752. M. J. Falk, B. R. Landry, and J. E. Subotnik, J. Phys. Chem. B 118, 8108 (2014). https://doi.org/10.1021/jp5011346
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–2422. T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100, 7884 (1996). https://doi.org/10.1021/jp953105a23. T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. A 101, 6389 (1997). https://doi.org/10.1021/jp970842t24. T. J. Martínez and R. D. Levine, J. Chem. Soc., Faraday Trans. 93, 941 (1997). https://doi.org/10.1039/a605958i 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. 5353. B. Mignolet and B. F. E. Curchod, J. Chem. Phys. 148, 134110 (2018). https://doi.org/10.1063/1.5022877 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-the-fly 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,5453. B. Mignolet and B. F. E. Curchod, J. Chem. Phys. 148, 134110 (2018). https://doi.org/10.1063/1.502287754. B. Mignolet and B. F. E. Curchod, J. Phys. Chem. A 123, 3582 (2019). https://doi.org/10.1021/acs.jpca.9b00940
3. 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 (NTBF), NTBF × (NTBF + 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.2626. B. F. E. Curchod, W. J. Glover, and T. J. Martínez, J. Phys. Chem. A 124, 6133 (2020). https://doi.org/10.1021/acs.jpca.0c04113 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 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. 2626. B. F. E. Curchod, W. J. Glover, and T. J. Martínez, J. Phys. Chem. A 124, 6133 (2020). https://doi.org/10.1021/acs.jpca.0c04113 (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.
4. 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 (Ninit) with the number of runs per initial condition (nrun), $NESTSH=nrun×Ninit$. 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 $NESAIMS(t)=∑k=1NinitNTBFk(t)×NTBFk(t)+1/2$, while for (E/O)SSAIMS it 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, $NES(E/O)SSAIMS(t)=∑j=1nrun∑k=1NinitNTBFk,j(t)×NTBFk,j(t)+1/2$. 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. 2525. B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martínez, Chem. Phys. 347, 3 (2008). https://doi.org/10.1016/j.chemphys.2008.01.014. 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.
B. Computational details
1. 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 self-consistent field (SA-CASSCF)55,5655. H.-J. Werner, Advances in Chemical Physics: Ab Initio Methods in Quantum Chemistry Part 2 (Wiley, 1987), Vol. 69, p. 1.56. B. O. Roos, Advances in Chemical Physics: Ab Initio Methods in Quantum Chemistry Part 2 (Wiley, 1987), Vol. 69, p. 399. and a 6-31G* basis set.57,5857. R. Ditchfield, W. J. Hehre, and J. A. Pople, J. Chem. Phys. 54, 724 (1971). https://doi.org/10.1063/1.167490258. P. C. Hariharan and J. A. Pople, Theor. Chem. Acc. 28, 213 (1973). https://doi.org/10.1007/bf00533485 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.5959. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, and M. Wang, Molpro, version 2012.1, a package of ab initio programs, 2012. 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 package5959. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, and M. Wang, Molpro, version 2012.1, a package of ab initio programs, 2012. (for the TSH dynamics) and TeraChem60–6460. I. S. Ufimtsev and T. J. Martínez, J. Chem. Theory Comput. 4, 222 (2008). https://doi.org/10.1021/ct700268q61. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput. 5, 3138 (2009). https://doi.org/10.1021/ct900433g62. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput. 5, 2619 (2009). https://doi.org/10.1021/ct900300463. S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, S. I. L. Kokkila-Schumacher, N. Luehr, J. W. Snyder, Jr., C. Song, A. V. Titov, I. S. Ufimtsev et al., J. Chem. Phys. 152, 224110 (2020). https://doi.org/10.1063/5.000761564. S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, C. M. Isborn, S. I. Kokkila-Schumacher, X. Li, F. Liu, N. Luehr, J. W. Snyder, Jr. et al., Wiley Interdiscip. Rev.: Comput. Mol. Sci. 11, e1494 (2020). https://doi.org/10.1002/wcms.1494 for the AIMS and (E/O)SSAIMS dynamics.
2. Nuclear dynamics
The AIMS and (E/O)SSAIMS dynamics for cyclopropanone and fulvene were performed with the AIMS/MOLPRO interface.2525. B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martínez, Chem. Phys. 347, 3 (2008). https://doi.org/10.1016/j.chemphys.2008.01.014 For 1,2-dithiane, the AIMS implementation in TeraChem was used.65,6665. J. W. Snyder, Jr., E. G. Hohenstein, N. Luehr, and T. J. Martínez, J. Chem. Phys. 143, 154107 (2015). https://doi.org/10.1063/1.493261366. J. W. Snyder, Jr., B. F. E. Curchod, and T. J. Martínez, J. Phys. Chem. Lett. 7, 2444 (2016). https://doi.org/10.1021/acs.jpclett.6b00970 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–6967. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, and L. González, J. Chem. Theory Comput. 7, 1253 (2011). https://doi.org/10.1021/ct100739468. S. Mai, P. Marquetand, and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 8, e1370 (2018). https://doi.org/10.1002/wcms.137069. S. Mai, M. Richter, M. Heindl, M. F. S. J. Menger, A. Atkins, M. Ruckenbauer, F. Plasser, M. Oppel, P. Marquetand, and L. González, SHARC2.0: Surface hopping including arbitrary couplings—Program package for non-adiabatic dynamics, sharc-md.org, 2018. 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 S1 population decay of ESSAIMS and (decoherence-corrected) TSH was in agreement: This criterion was fulfilled with five runs for cyclopropanone, seven for fulvene, and eight for 1,2-dithiane (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.7070. F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand, and L. González, J. Chem. Theory Comput. 12, 1207 (2016). https://doi.org/10.1021/acs.jctc.5b01148 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,7149. G. Granucci, M. Persico, and A. Zoccante, J. Chem. Phys. 133, 134111 (2010). https://doi.org/10.1063/1.348900471. G. Granucci and M. Persico, J. Chem. Phys. 126, 134114 (2007). https://doi.org/10.1063/1.2715585 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. 4949. G. Granucci, M. Persico, and A. Zoccante, J. Chem. Phys. 133, 134111 (2010). https://doi.org/10.1063/1.3489004. 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 in the following). We note that the original implementation of the EDC in SHARC was employed, where the amplitudes of the non-running states are damped, instead of the populations, as done, for example, in Newton-X.49,72,7349. G. Granucci, M. Persico, and A. Zoccante, J. Chem. Phys. 133, 134111 (2010). https://doi.org/10.1063/1.348900472. M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico, and H. Lischka, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 4, 26 (2014). https://doi.org/10.1002/wcms.115873. A. Prlj, NEWTON-X forum: Implementation of decoherence correction, 2019. All the dTSH trajectories strictly satisfy the internal consistency criterion discussed in Ref. 7171. G. Granucci and M. Persico, J. Chem. Phys. 126, 134114 (2007). https://doi.org/10.1063/1.2715585, 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.
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–3027. G. Cui, Y. Ai, and W. Fang, J. Phys. Chem. A 114, 730 (2010). https://doi.org/10.1021/jp908936u28. G. Cui and W. Fang, J. Phys. Chem. A 115, 1547 (2011). https://doi.org/10.1021/jp110632g29. M. Filatov, S. K. Min, and C. H. Choi, Phys. Chem. Chem. Phys. 21, 2489 (2019). https://doi.org/10.1039/c8cp07104g30. J. Suchan, J. Janoš, and P. Slavíček, J. Chem. Theory Comput. 16, 5809 (2020). https://doi.org/10.1021/acs.jctc.0c00512 All these previous studies highlight the fact that the S1 population decay of cyclopropanone following photoexcitation is rather straightforward—the nuclear wavepacket decays from S1 to S0 without any significant back population transfer. The S1 decay has been characterized by a latency time in S1 of ∼25 fs, followed by a stepwise decay of the S1 population.3030. J. Suchan, J. Janoš, and P. Slavíček, J. Chem. Theory Comput. 16, 5809 (2020). https://doi.org/10.1021/acs.jctc.0c00512 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–2927. G. Cui, Y. Ai, and W. Fang, J. Phys. Chem. A 114, 730 (2010). https://doi.org/10.1021/jp908936u28. G. Cui and W. Fang, J. Phys. Chem. A 115, 1547 (2011). https://doi.org/10.1021/jp110632g29. M. Filatov, S. K. Min, and C. H. Choi, Phys. Chem. Chem. Phys. 21, 2489 (2019). https://doi.org/10.1039/c8cp07104g 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.
1. AIMS vs (E/O)SSAIMS
We start by comparing the S1 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 S1 population starts decaying after 25 fs, before plateauing at 70% of S1 population. ESSAIMS (ϵ = 10−10 a.u.) remains stable enough to describe the beginning of the subsequent S1 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 S1 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 S0 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 S1 population decay (see the blue curve in Fig. 3). The population trace obtained with this threshold remains within the error bars of AIMS and ESSAIMS (ϵ = 10−10 a.u.). The average number of TBFs does not surpass 1.25 for the ESSAIMS (ϵ = 10−5 a.u.) dynamics—well under the average number of ESSAIMS (ϵ = 10−10 a.u.)—and stays below 1.1 for the largest part of the dynamics. As observed earlier, the average number of TBFs rises initially but drops to almost 1.0 when the first plateau is reached. It then increases again every time new TBFs are spawned, permitting the population transfer to S0.
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 S1 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.
2. 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 S1 population decay for the first 75 fs of dynamics. The population decay starts after 25 fs and plateaus at around 65% of S1 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 S1 is not dampened, artificially increasing the probabilities of hops back to the S1 state, leading to an overall slowdown of the S1 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 S1 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.3535. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f
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 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 S1, fulvene can encounter two different conical intersections:31–3531. M. J. Bearpark, F. Bernardi, M. Olivucci, M. A. Robb, and B. R. Smith, J. Am. Chem. Soc. 118, 5254 (1996). https://doi.org/10.1021/ja954279932. S. Alfalah, S. Belz, O. Deeb, M. Leibscher, J. Manz, and S. Zilberg, J. Chem. Phys. 130, 124318 (2009). https://doi.org/10.1063/1.308954633. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. J. Bearpark, and M. A. Robb, Phys. Chem. Chem. Phys. 12, 15725 (2010). https://doi.org/10.1039/c0cp01757d34. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. A. Robb, and M. J. Bearpark, J. Chem. Phys. 137, 22A548 (2012). https://doi.org/10.1063/1.476508735. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f a peaked conical intersection leading to an efficient deactivation to S0 mediated by a twist of the C=CH2 moiety or a strongly sloped conical intersection reached by a stretching of the C=CH2 bond, leading to a transfer of the molecule to S0 but also a possible reflection process bringing the molecule back to the S1 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. 3535. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f for more information).
The bottom panel of Fig. 5 shows the S1 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 S1 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 S1 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. 3535. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f 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.
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:3636. C. D. Rankine, J. P. F. Nunes, M. S. Robinson, P. D. Lane, and D. A. Wann, Phys. Chem. Chem. Phys. 18, 27170 (2016). https://doi.org/10.1039/c6cp05518d Upon photoexcitation, dithiane commences an ultrafast ring-opening process in S1 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 work7474. A. B. Stephansen, R. Y. Brogaard, T. S. Kuhlman, L. B. Klein, J. B. Christensen, and T. I. Sølling, J. Am. Chem. Soc. 134, 20279 (2012). https://doi.org/10.1021/ja310540a). 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 S1 population begins to decay after 25 fs (bottom panel in Fig. 6), leading to an S1 population drop to 40% within 60 fs. Subsequently, the S1 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 S1 population is overestimated during the revival process. After 100 fs, the S1 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 S1 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 S1 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 S1 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 S1 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 S1 population decay, monitoring the S0 population reveals larger deviations between the methods (Fig. 7). OSSAIMS and AIMS show that some population appears in S0 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 S0 population takes place earlier (mirroring the behavior observed for the S1 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 S0 population dynamics up to 100 fs of dynamics (light blue dashed curve in Fig. 7). However, the S0 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.7575. We note that a threshold of ϵ = 10−4 a.u. does not improve the dynamics significantly (see supplementary material).
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 S0 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 S0 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,2-dithiane 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 off-diagonal 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.
In summary, we applied the novel framework of stochastic-selection 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 S1 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 S1/S0 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 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.
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.
This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 803718, project SINDAM). L.M.I. acknowledges the EPSRC for an EPSRC Doctoral Studentship (Grant No. EP/R513039/1). T.J.M. acknowledges support from the Chemical Sciences, Geosciences and Biosciences Division of the Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy. This work made use of the facilities of the Hamilton HPC Service of Durham University and benefited from workshops organized by the E-CAM [European Union’s Horizon 2020 research and innovation program (Grant No. 676531)].
The data that support the findings of this study are available within the article and its supplementary material.
1. 1. J. C. Tully, Theor. Chem. Acc. 103, 173 (2000). https://doi.org/10.1007/978-3-662-10421-7_3, Google ScholarCrossref
2. 2. M. Born and R. Oppenheimer, Ann. Phys. 389, 457 (1927). https://doi.org/10.1002/andp.19273892002, Google ScholarCrossref
3. 3. H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990). https://doi.org/10.1016/0009-2614(90)87014-i, Google ScholarCrossref
4. 4. M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, Phys. Rep. 324, 1 (2000). https://doi.org/10.1016/S0370-1573(99)00047-2, Google ScholarCrossref
5. 5. G. A. Worth, H. D. Meyer, and L. S. Cederbaum, in Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, Advanced Series in Physical Chemistry Vol. 15, edited by W. Domcke, D. R. Yarkony, and H. Köppel (World Scientific, 2004), Chap. 14, pp. 583–617. Google ScholarCrossref
6. 6. H.-D. Meyer, F. Gatti, and G. A. Worth, Multidimensional Quantum Dynamics (John Wiley & Sons, 2009). Google ScholarCrossref
7. 7. G. W. Richings and S. Habershon, J. Chem. Phys. 148, 134116 (2018). https://doi.org/10.1063/1.5024869, Google ScholarScitation, ISI
8. 8. G. W. Richings and S. Habershon, J. Phys. Chem. A 124, 9299 (2020). https://doi.org/10.1021/acs.jpca.0c06125, Google ScholarCrossref
9. 9. M. Barbatti, Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 620 (2011). https://doi.org/10.1002/wcms.64, Google ScholarCrossref
10. 10. J. C. Tully, J. Chem. Phys. 137, 22A301 (2012). https://doi.org/10.1063/1.4757762, Google ScholarScitation, ISI
11. 11. J. C. Tully, J. Chem. Phys. 93, 1061 (1990). https://doi.org/10.1063/1.459170, Google ScholarScitation, ISI
12. 12. S. Hammes-Schiffer and J. C. Tully, J. Chem. Phys. 101, 4657 (1994). https://doi.org/10.1063/1.467455, Google ScholarScitation, ISI
13. 13. R. Crespo-Otero and M. Barbatti, Chem. Rev. 118, 7026 (2018). https://doi.org/10.1021/acs.chemrev.7b00577, Google ScholarCrossref
14. 14. F. Agostini, S. K. Min, A. Abedi, and E. K. U. Gross, J. Chem. Theory Comput. 12, 2127 (2016). https://doi.org/10.1021/acs.jctc.5b01180, Google ScholarCrossref
15. 15. S. K. Min, F. Agostini, I. Tavernelli, and E. K. U. Gross, J. Phys. Chem. Lett. 8, 3048 (2017). https://doi.org/10.1021/acs.jpclett.7b01249, Google ScholarCrossref
16. 16. G. A. Worth, M. A. Robb, and I. Burghardt, Faraday Discuss. 127, 307 (2004). https://doi.org/10.1039/b314253a, Google ScholarCrossref
17. 17. B. Lasorne, M. J. Bearpark, M. A. Robb, and G. A. Worth, Chem. Phys. Lett. 432, 604 (2006). https://doi.org/10.1016/j.cplett.2006.10.099, Google ScholarCrossref
18. 18. G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, Int. Rev. Phys. Chem. 34, 269 (2015). https://doi.org/10.1080/0144235x.2015.1051354, Google ScholarCrossref
19. 19. D. V. Shalashilin, J. Chem. Phys. 130, 244101 (2009). https://doi.org/10.1063/1.3153302, Google ScholarScitation, ISI
20. 20. K. Saita and D. V. Shalashilin, J. Chem. Phys. 137, 22A506 (2012). https://doi.org/10.1063/1.4734313, Google ScholarScitation, ISI
21. 21. D. V. Makhov, C. Symonds, S. Fernandez-Alberti, and D. V. Shalashilin, Chem. Phys. 493, 200 (2017). https://doi.org/10.1016/j.chemphys.2017.04.003, Google ScholarCrossref
22. 22. T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100, 7884 (1996). https://doi.org/10.1021/jp953105a, Google ScholarCrossref
23. 23. T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. A 101, 6389 (1997). https://doi.org/10.1021/jp970842t, Google ScholarCrossref
24. 24. T. J. Martínez and R. D. Levine, J. Chem. Soc., Faraday Trans. 93, 941 (1997). https://doi.org/10.1039/a605958i, Google ScholarCrossref
25. 25. B. G. Levine, J. D. Coe, A. M. Virshup, and T. J. Martínez, Chem. Phys. 347, 3 (2008). https://doi.org/10.1016/j.chemphys.2008.01.014, Google ScholarCrossref
26. 26. B. F. E. Curchod, W. J. Glover, and T. J. Martínez, J. Phys. Chem. A 124, 6133 (2020). https://doi.org/10.1021/acs.jpca.0c04113, Google ScholarCrossref
27. 27. G. Cui, Y. Ai, and W. Fang, J. Phys. Chem. A 114, 730 (2010). https://doi.org/10.1021/jp908936u, Google ScholarCrossref
28. 28. G. Cui and W. Fang, J. Phys. Chem. A 115, 1547 (2011). https://doi.org/10.1021/jp110632g, Google ScholarCrossref
29. 29. M. Filatov, S. K. Min, and C. H. Choi, Phys. Chem. Chem. Phys. 21, 2489 (2019). https://doi.org/10.1039/c8cp07104g, Google ScholarCrossref
30. 30. J. Suchan, J. Janoš, and P. Slavíček, J. Chem. Theory Comput. 16, 5809 (2020). https://doi.org/10.1021/acs.jctc.0c00512, Google ScholarCrossref
31. 31. M. J. Bearpark, F. Bernardi, M. Olivucci, M. A. Robb, and B. R. Smith, J. Am. Chem. Soc. 118, 5254 (1996). https://doi.org/10.1021/ja9542799, Google ScholarCrossref
32. 32. S. Alfalah, S. Belz, O. Deeb, M. Leibscher, J. Manz, and S. Zilberg, J. Chem. Phys. 130, 124318 (2009). https://doi.org/10.1063/1.3089546, Google ScholarScitation, ISI
33. 33. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. J. Bearpark, and M. A. Robb, Phys. Chem. Chem. Phys. 12, 15725 (2010). https://doi.org/10.1039/c0cp01757d, Google ScholarCrossref
34. 34. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. A. Robb, and M. J. Bearpark, J. Chem. Phys. 137, 22A548 (2012). https://doi.org/10.1063/1.4765087, Google ScholarScitation, ISI
35. 35. L. M. Ibele and B. F. E. Curchod, Phys. Chem. Chem. Phys. 22, 15183 (2020). https://doi.org/10.1039/d0cp01353f, Google ScholarCrossref
36. 36. C. D. Rankine, J. P. F. Nunes, M. S. Robinson, P. D. Lane, and D. A. Wann, Phys. Chem. Chem. Phys. 18, 27170 (2016). https://doi.org/10.1039/c6cp05518d, Google ScholarCrossref
37. 37. B. F. E. Curchod and T. J. Martínez, Chem. Rev. 118, 3305 (2018). https://doi.org/10.1021/acs.chemrev.7b00423, Google ScholarCrossref
38. 38. L. M. Ibele, A. Nicolson, and B. F. E. Curchod, Mol. Phys. 118, e1665199 (2020). https://doi.org/10.1080/00268976.2019.1665199, Google ScholarCrossref
39. 39. J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016). https://doi.org/10.1146/annurev-physchem-040215-112245, Google ScholarCrossref
40. 40. F. Agostini and B. F. E. Curchod, Wiley Interdiscip. Rev. Comput. Mol. Sci. 9, e1417 (2019). https://doi.org/10.1002/wcms.1417, Google ScholarCrossref
41. 41. E. R. Bittner and P. J. Rossky, J. Chem. Phys. 103, 8130 (1995). https://doi.org/10.1063/1.470177, Google ScholarScitation, ISI
42. 42. M. Thachuk, M. Y. Ivanov, and D. M. Wardlaw, J. Chem. Phys. 109, 5747 (1998). https://doi.org/10.1063/1.477197, Google ScholarScitation, ISI
43. 43. J.-Y. Fang and S. Hammes-Schiffer, J. Phys. Chem. A 103, 9399 (1999). https://doi.org/10.1021/jp991602b, Google ScholarCrossref
44. 44. A. W. Jasper, S. Nangia, C. Zhu, and D. G. Truhlar, Acc. Chem. Res. 39, 101 (2006). https://doi.org/10.1021/ar040206v, Google ScholarCrossref
45. 45. J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 024105 (2011). https://doi.org/10.1063/1.3506779, Google ScholarScitation, ISI
46. 46. J. E. Subotnik and N. Shenvi, J. Chem. Phys. 134, 244114 (2011). https://doi.org/10.1063/1.3603448, Google ScholarScitation, ISI
47. 47. M. Persico and G. Granucci, Theor. Chem. Acc. 133, 1526 (2014). https://doi.org/10.1007/s00214-014-1526-1, Google ScholarCrossref
48. 48. B. F. E. Curchod and I. Tavernelli, J. Chem. Phys. 138, 184112 (2013). https://doi.org/10.1063/1.4803835, Google ScholarScitation, ISI
49. 49. G. Granucci, M. Persico, and A. Zoccante, J. Chem. Phys. 133, 134111 (2010). https://doi.org/10.1063/1.3489004, Google ScholarScitation, ISI
50. 50. H. M. Jaeger, S. Fischer, and O. V. Prezhdo, J. Chem. Phys. 137, 22A545 (2012). https://doi.org/10.1063/1.4757100, Google ScholarScitation, ISI
51. 51. N. Shenvi and W. Yang, J. Chem. Phys. 137, 22A528 (2012). https://doi.org/10.1063/1.4746407, Google ScholarScitation, ISI
52. 52. M. J. Falk, B. R. Landry, and J. E. Subotnik, J. Phys. Chem. B 118, 8108 (2014). https://doi.org/10.1021/jp5011346, Google ScholarCrossref
53. 53. B. Mignolet and B. F. E. Curchod, J. Chem. Phys. 148, 134110 (2018). https://doi.org/10.1063/1.5022877, Google ScholarScitation, ISI
54. 54. B. Mignolet and B. F. E. Curchod, J. Phys. Chem. A 123, 3582 (2019). https://doi.org/10.1021/acs.jpca.9b00940, Google ScholarCrossref
55. 55. H.-J. Werner, Advances in Chemical Physics: Ab Initio Methods in Quantum Chemistry Part 2 (Wiley, 1987), Vol. 69, p. 1. Google Scholar
56. 56. B. O. Roos, Advances in Chemical Physics: Ab Initio Methods in Quantum Chemistry Part 2 (Wiley, 1987), Vol. 69, p. 399. Google ScholarCrossref
57. 57. R. Ditchfield, W. J. Hehre, and J. A. Pople, J. Chem. Phys. 54, 724 (1971). https://doi.org/10.1063/1.1674902, Google ScholarScitation, ISI
58. 58. P. C. Hariharan and J. A. Pople, Theor. Chem. Acc. 28, 213 (1973). https://doi.org/10.1007/bf00533485, Google ScholarCrossref
59. 59. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, and M. Wang, Molpro, version 2012.1, a package of ab initio programs, 2012. Google Scholar
60. 60. I. S. Ufimtsev and T. J. Martínez, J. Chem. Theory Comput. 4, 222 (2008). https://doi.org/10.1021/ct700268q, Google ScholarCrossref
61. 61. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput. 5, 3138 (2009). https://doi.org/10.1021/ct900433g, Google ScholarCrossref
62. 62. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput. 5, 2619 (2009). https://doi.org/10.1021/ct9003004, Google ScholarCrossref
63. 63. S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, S. I. L. Kokkila-Schumacher, N. Luehr, J. W. Snyder, Jr., C. Song, A. V. Titov, I. S. Ufimtsev et al., J. Chem. Phys. 152, 224110 (2020). https://doi.org/10.1063/5.0007615, Google ScholarScitation, ISI
64. 64. S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, C. M. Isborn, S. I. Kokkila-Schumacher, X. Li, F. Liu, N. Luehr, J. W. Snyder, Jr. et al., Wiley Interdiscip. Rev.: Comput. Mol. Sci. 11, e1494 (2020). https://doi.org/10.1002/wcms.1494, Google ScholarCrossref
65. 65. J. W. Snyder, Jr., E. G. Hohenstein, N. Luehr, and T. J. Martínez, J. Chem. Phys. 143, 154107 (2015). https://doi.org/10.1063/1.4932613, Google ScholarScitation, ISI
66. 66. J. W. Snyder, Jr., B. F. E. Curchod, and T. J. Martínez, J. Phys. Chem. Lett. 7, 2444 (2016). https://doi.org/10.1021/acs.jpclett.6b00970, Google ScholarCrossref
67. 67. M. Richter, P. Marquetand, J. González-Vázquez, I. Sola, and L. González, J. Chem. Theory Comput. 7, 1253 (2011). https://doi.org/10.1021/ct1007394, Google ScholarCrossref
68. 68. S. Mai, P. Marquetand, and L. González, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 8, e1370 (2018). https://doi.org/10.1002/wcms.1370, Google ScholarCrossref
69. 69. S. Mai, M. Richter, M. Heindl, M. F. S. J. Menger, A. Atkins, M. Ruckenbauer, F. Plasser, M. Oppel, P. Marquetand, and L. González, SHARC2.0: Surface hopping including arbitrary couplings—Program package for non-adiabatic dynamics, sharc-md.org, 2018. Google Scholar
70. 70. F. Plasser, M. Ruckenbauer, S. Mai, M. Oppel, P. Marquetand, and L. González, J. Chem. Theory Comput. 12, 1207 (2016). https://doi.org/10.1021/acs.jctc.5b01148, Google ScholarCrossref
71. 71. G. Granucci and M. Persico, J. Chem. Phys. 126, 134114 (2007). https://doi.org/10.1063/1.2715585, Google ScholarScitation, ISI
72. 72. M. Barbatti, M. Ruckenbauer, F. Plasser, J. Pittner, G. Granucci, M. Persico, and H. Lischka, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 4, 26 (2014). https://doi.org/10.1002/wcms.1158, Google ScholarCrossref
73. 73. A. Prlj, NEWTON-X forum: Implementation of decoherence correction, 2019. Google Scholar
74. 74. A. B. Stephansen, R. Y. Brogaard, T. S. Kuhlman, L. B. Klein, J. B. Christensen, and T. I. Sølling, J. Am. Chem. Soc. 134, 20279 (2012). https://doi.org/10.1021/ja310540a, Google ScholarCrossref
75. 75.We note that a threshold of ϵ = 10−4 a.u. does not improve the dynamics significantly (see supplementary material).