Critical Assessment of Curvature-Driven Surface Hopping Algorithms

Critical Assessment of Curvature-Driven Surface Hopping Algorithms TomásJíra, JiríJanos, and Petr Slavícek* Cite This: J. Chem. Theory Comput. 2025, 21, 9784−9798 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: Trajectory surface-hopping (TSH) methods have become the most used approach in nonadiabatic molecular dynamics. The increasingly popular curvature-driven schemes represent a subset of TSH based on the implicit local diabatization of potential energy surfaces. Their appeal partly stems from compatibility with machinelearning frameworks that often provide only local PES information. Here, we critically assess the limitations of these curvature-based algorithms by examining three challenging scenarios: (i) dynamics involving more than two strongly coupled electronic states; (ii) trivial crossings; and (iii) spurious transitions arising from small discontinuities in multireference potential energy surfaces. Furthermore, we extend the Landau−Zener surface hopping (LZSH) method beyond two-state systems and introduce practical modifications to enhance its robustness. The performance is benchmarked on both lowand higher-dimensional model Hamiltonians, as well as realistic molecular systems treated with ab initio methods. While curvaturedriven TSH using the explicit electronic coefficient propagation qualitatively captures the dynamics in most cases, we find no regime where it outperforms LZSH, especially when trivial crossings, multistate crossings, or discontinuities are encountered. Hence, we advocate for using a conceptually simple but solid LZSH method when nonadiabatic couplings are not available. 1. INTRODUCTION Molecular dynamics has proven to be a powerful tool for quantitatively modeling chemical reactions. 1 For photochemical systems, however, simulations must typically go beyond the adiabatic approximation and account for coupled electronic and nuclear motion. 2 The dynamics in such cases appear on multiple potential energy surfaces (PESs), with coupling treated beyond perturbation theory. With suitably chosen model potentials, highly accurate quantum approaches such as the multiconfigurational timedependent Hartree (MCTDH) method can be employed. 3,4 More often, however, we study reactions without detailed prior knowledge, relying instead on on-the-fly electronic-structure data. 5,6 As these data are local, trajectory-based methods incorporating surface hopping between electronic states have become standard for simulating nonadiabatic dynamics, beginning with the pioneering ab initio studies. 2,7,8 These methods integrate seamlessly with ab initio codes, offer intuitive interpretations of dynamics, and are computationally efficient and trivially parallelizable. A key distinction between surface hopping schemes lies in how electronic transitions are modeled. For many years, the field was dominated by the Landau−Zener−Stuckelberg− Majorana framework (here simply referred to as Landau− Zener, or LZ). 9−13 This formalism assumes two-state crossings with linear (diabatic) potentials, constant coupling and constant velocity, yielding reliable predictions when only two states are coupled. It has proven valuable in both qualitative and quantitative contexts, particularly in charge-transfer reactions. 14 Notably, the LZ model was also employed in the first molecular dynamics surface hopping treatments. 15 To address the limitations of the LZ approach in multistate systems, Tully introduced the fewest-switches surface hopping (FSSH) algorithm. 16 Designed for general potential surfaces without localizing the crossing region, FSSH gained popularity alongside advances in electronic structure methods and quickly became a standard in nonadiabatic dynamics. 17 Yet several issues persist. Chief among them is the lack of electronic decoherence in FSSH, typically handled through ad hoc corrections, 18,19 though more principled approaches have emerged. 20−23 FSSH also struggles with trivial crossings (degenerate diabatic states with negligible coupling) common e.g., in molecular crystals. 24−27 Moreover, the algorithm is often implemented�due to its simplicity�using nonadiabatic coupling vectors (NACVs), which are often divergent or unavailable. As an alternative, some implementations directly Received: July 17, 2025 Revised: August 28, 2025 Accepted: August 28, 2025 Published: September 16, 2025 Articlepubs.acs.org/JCTC © 2025 The Authors. Published by American Chemical Society 9784 https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 This article is licensed under CC-BY 4.0 Downloaded via UNIV CHEMISTRY & TECHNOLOGY PRAGUE on October 23, 2025 at 14:06:21 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. approximate the time derivative coupling (TDC) via electronic wave function overlaps, 28,29 enhancing FSSH stability. These limitations have renewed interest in alternative, often more pragmatic, hopping schemes. 30−36 This trend reflects a growing consensus that the choice of electronic structure method more strongly determines simulation outcomes than the hopping algorithm itself. 37−44 Consequently, robustness and computational feasibility are increasingly prioritized. While simple schemes that switch states upon encountering a new surface are sometimes used, 15 more physical approaches revisit the two-state problem. For instance, the mapping approach to surface hopping (MASH) addresses decoherence rigorously for two-state systems, 20 and Landau−Zener surface hopping (LZSH) has been reformulated in the adiabatic representation. 45 The latter removes the need for NACVs and allows a natural treatment of intersystem crossings, which is used in numerous studies. 31−33,46,47 More recently, Baeck and An proposed an approximation to nonadiabatic couplings based on PES curvature, 48 leading to curvature-driven fewestswitches schemes that allow transitions throughout a trajectory. 30,40,49−51 These methods can be easily combined with any electronic structure code capable of delivering potential energy surfaces, without the need to explicitly calculate the nonadiabatic couplings. The present work focuses on curvature-based approaches derived under the two-state assumption. This includes both the LZSH and curvature-driven fewest switches surface hopping methods. We aim to identify regimes where such methods fail in photodynamics simulations. We consider three critical challenges: (i) three-state interactions, where methods designed for two states may misrepresent couplings, (ii) trivial crossings, which occur frequently in condensed phase systems and are known to cause failures in conventional hopping methods, 52,53 and (iii) slight discontinuities in PESs, common in ab initio surfaces, which may degrade curvature-based schemes. The remainder of this article is organized as follows. Section 2introduces the theoretical framework and algorithms. Section 3provides computational details for all the calculations. Section 4.1 presents benchmark results on model potentials designed to stress-test the methods. Section 3.2 showcases the surface hopping algorithms on a well-known vibronic coupling model of the uracil cation. Sections 4.3 and 4.4 extend this analysis to realistic molecular systems, focusing on the effects of PES continuity. 2. THEORY 2.1. Trajectory Surface Hopping. In this section, we outline the theoretical framework of trajectory surface hopping (TSH). The key idea underlying TSH is that the nuclei are treated classically and propagated via Newton’s equations of motion, whereas the electrons are treated quantum mechanically. Under this approximation, the nuclei evolve on a potential energy surface obtained using the electronic timeindependent Schrodinger equation. The electronic wave function Φ(r,t;R(t)) is then expanded as a linear combination of adiabatic eigenstates ϕi(r;R(t)) at nuclear position R(t) of the electronic Hamiltonian as =t t c t r tr R R( , ; ( )) ( ) ( ; ( )) k kk (1) where the time-dependent coefficients ck(t) carry information about the population and phase of each electronic state ϕk. For simplicity, we will omit the explicit dependence of the wave functions on their dynamical variables in subsequent equations. To obtain the equations of motion for the coefficients ck, we substitute eq 1 into the time-dependent Schrodinger equation (TDSE), obtaining =ic tc H i t d d j k k jk j kel i k j j j j j j y { z z z z z z (2) where the term = jk jt k is the time derivative coupling. In the adiabatic basis, the electronic Hamiltonian Hel is diagonal, so the TDC alone governs the transfer of population between adiabatic states. Because quantum chemical calculations can provide NACVs, =djk jR k , instead of time derivative couplings, we can employ the chain rule to write = = · = · t t R R v d d d jk j k j k jk (3) With this manipulation, we propagate the nuclei together with the electronic wave function on a single adiabatic surface, relying only on electronic gradients and NACVs. The TSH algorithm needs to be supplemented by a procedure to allow transitions (hops) between different (adiabatic) surfaces. Several methods are used, including a simple approach in which hops occur whenever the energy gap between adiabatic surfaces appears below a preset threshold. 54 Methods used within this document are described in the following sections. 2.2. Landau−Zener Surface Hopping. The Landau− Zener surface hopping algorithm, originally developed nearly a century ago, provides a remarkably straightforward and powerful framework for modeling nonadiabatic transitions. 9−11 In the LZSH framework, an analytic solution of TDSE in eq 2 for a simplified model is applied throughout the dynamics. Hence, the electronic TDSE is not propagated along trajectories. Consider a two-level system described by the diabatic Hamiltonian representing the crossing between potential energy surfaces = t H H t H 12 21 Ä Ç Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ (4) where αis an arbitrary constant and the diabatic coupling H12 is independent of time. The corresponding TDSE (eq2in diabatic basis) is written as =i t c t c t c t c t H d d ( ) ( ) ( ) ( ) 1 2 1 2 Ä Ç Å Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ä Ç Å Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ (5) with c1and c2representing the coefficients of the diabatic wave function. Under the initial condition c1(t→ −∞) = 1, an exact solution of eq 5 leads to =c t H ( ) exp 112 2 i k j j j j j y { z z z z z (6) The hopping probability from state 1 to state 2 then reads 55,56 Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9785 = = = P c t c t H ( ) 1 ( ) 1 exp 2 1 2 LZ,dia 2212 12 2 i k j j j j j y { z z z z z (7) After simple algebraic manipulation and switching to general electronic state indexes, the adiabatic transition probability is obtained as 45 =P Z Z exp 2 i k ik ik LZ 3 i k j j j j j j j y { z z z z z z z (8) where Zik =|Ek−Ei|is the energy difference between states with Eiand Ekbeing the energies of adiabatic electronic states i and jat the current nuclear geometry. This probability is evaluated upon encountering the minimum of the energy difference. It is important to note that the LZSH algorithm is intrinsically limited by its assumptions. It is formulated solely for two-state systems with linearly time-dependent energies and constant diabatic coupling. Although this constitutes a considerable approximation, the method provides a basis for our qualitative understanding of nonadiabatic processes and often yields quantitative results that are in reasonable agreement with those obtained using the FSSH approach. 32,37,57 Since the LZSH algorithm is defined only for two-state crossings, we have extended the algorithm for multistate systems. When a multistate crossing is identified, our modified LZSH algorithm computes the transition probabilities from the current state only to a state for which Zik has a maximum value. The motivation behind the scheme stems from the fact that the two interacting states will exhibit the largest second derivatives, while the noninteracting states are expected to have relatively smaller second derivatives. We will refer to this algorithm as LZSHmaxcurv. We also inspect a variant of LZSH, where we only calculate the transition probability to the energetically closest state. This modification will be referred to as LZSHnearest. The two versions of multistate LZSH are schematically depicted in Figure 1. Moreover, we investigated a version where we allow hops only to the state with maximum hopping probability, but it provides the same results as the LZSHnearest algorithm. Hence, we do not discuss it further. Note that for a two-state system, all of these multistate variants reduce to the original LZSH scheme. An alternative multistate extension of the Landau−Zener model has been proposed by Smith and Akimov, in which no final state is favored but the hopping probabilities are rather normalized accordingly. 58 We deliberately chose not to employ this variant, as it yields incorrect results for our specific threestate model systems. Smith and Akimov further suggested interpolating the adiabatic energy curves to improve the accuracy of the second derivatives, which we also do not use. 58 An additional weakness of LZSH, and generally of algorithms based on PES’s curvature, stems from encountering discontinuities in electronic energies along trajectories. Such discontinuities are ubiquitous in multireference electronic structure methods and are caused by either orbital rotations into the active space or state flipping in the state-averaging procedure. While such discontinuities should be ideally avoided, general consensus tolerates “reasonable” discontinuities, providing that the character of the active electronic state has not changed. Such discontinuities then create spurious (and often very sharp) minima in Zik. The LZSH algorithm spots such minima and, unable to recognize discontinuity from a true minimum, assigns hopping probabilities. Induced spurious hops then plague the dynamics. To mitigate the issue, we developed a simple scheme to detect spurious minima stemming from discontinuities in PESs. If LZSH finds a minimum in Zik, it compares numerical second Figure 1. Top row: Schematic trajectories for typical situations encountered in photodynamics depicting principles of the two multistate LZSH variants. While both schemes are identical when only the neighboring states are interacting, they differ if the interacting states encompass a noninteracting state. Bottom row: Illustration of the performance of the blocking algorithm for discontinuity-induced hops in three different situations encountered in typical trajectories. For typical avoided crossing (left), the αvalue is small, and the algorithm allows all hops. If a very sharp conical intersection is encountered (middle), values of αaround one are typical. In such a situation, the algorithm is not able to recognize from the known data whether a conical intersection or a discontinuity was encountered. Hence, a warning is issued. If the coefficient αexceeds 1.3 (situation on the right), a discontinuity almost certainly occurred, and the hop is blocked. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9786 derivatives of Zik at the minimum (Zik,min) and at the previous step (Zik,prev). In the case of an avoided crossing, these second derivatives should vary smoothly. However, if a discontinuity is encountered, an abrupt change in the second derivative is expected. Therefore, we define a ratio α=|(Zik,min −Zik,prev)/ Zik,min|. If α< 0.3, we consider PESs continuous. If α∈[0.3, 1.3], the code issues a warning to check the PESs manually. The numerical thresholds were set empirically. Note that if a sharp conical intersection is hit directly, values of α≈1 are typically observed. If α> 1.3, the curvature change is too large −a typical feature of a discontinuity. The blocking algorithm is visually depicted in Figure 1. Note that more advanced schemes can also be developed, e.g., if one more step further is taken, yet this requires additional computational effort. A comparison between simulations performed with and without the discontinuity correction for the LZSH method is provided in the Supporting Information. 2.3. Fewest Switches Surface Hopping. Let us now move beyond the simple Landau−Zener model and introduce the fewest switches surface hopping approach. Originally proposed by Tully, FSSH is among the most widely used algorithms for describing transitions between electronic states. 16 Its key objective is to minimize the total number of hops while ensuring that state populations remain consistent with the quantum mechanical probabilities, hence the term “fewest switches”. The hopping probability from state ito state kis given as =| | [ ] [ ]Pt c c c H c c t 2 i k i k i ki k i k iFS 2 i k j j j j j j y { z z z z z z (9) In the adiabatic basis, Hki = 0 for k≠iand the first term in the above equation vanishes. The transition probability, therefore, simplifies to =| | [ ]Pt c c c t 2 i k i k i k i FS 2 (10) If this probability evaluates to a negative value, it is set to zero. Upon a successful hop from ito k, the nuclear coordinates continue to evolve classically on the kth adiabatic potential energy surface, and the velocity is rescaled (typically along the direction of the nonadiabatic coupling vector) to conserve the total energy. In certain instances, the system attempts a hop but lacks the necessary kinetic energy to transition to the new potential surface. These occurrences are known as “frustrated” hops. While there are various ad hoc approaches to handling such events, we follow Tully’s original method, where the hops are simply disregarded, since the simulations in this work are set up so that the particles always have enough kinetic energy to jump. A known limitation of FSSH is that it does not inherently satisfy the internal consistency condition, implying that at any point in time, the fraction of trajectories on a given electronic state may not match the corresponding quantum probability. To address this issue, Granucci and Persico introduced an ad hoc decoherence correction that damps the electronic wave function coefficients, thereby restoring consistency. 18 We adopt this correction in our simulations, and the specific parameters for decoherence will be detailed in subsequent sections. 2.4. Baeck−An Scheme, κTDC, and κFSSH. A key limitation of the FSSH algorithm is its need for TDC, which are typically computed from NACVs −an approach that can be too expensive, unstable, or even impossible for some electronic structure methods. While overlaps between electronic wave functions offer an alternative, this method requires careful phase tracking and is not supported in most standard quantum chemistry codes. To address this challenge, Baeck and An proposed an approximation for the magnitude of the NACV in the vicinity of the crossing point that relies solely on the curvature of the potential energy surface 48 =dZ RZ 1 2 d d 1 ik ik ik BA 2 2 (11) where Zik =|Ek−Ei|. We can use the chain rule and express the TDC according to Baeck and An, denoted as σik λ, as = Z t R R Z t Z 1 2 d d d d 1 ik ik ik ik 2 2 i k j j j j j y { z z z z z (12) This expression is unsigned, due to the nature of the Baeck− An coupling. Should we now take the approximation that =0 Z t d d ik at the energy crossing (note that the Baeck−An coupling was derived for the energy crossing and its validity beyond it is only approximate 48 ), we arrive at the commonly used expression for TDC 30 = Z tZ 1 2 d d 1 ik ik ik 2 2 (13) with the convention that for k>ithe expression is positive, while for k<ione sets σik κ=−σki κ. If the argument of the square root in eq 13 or eq 12 is negative, the element of the TDC is set to zero. The expression in eq 13 is often referred to as κTDC. A FSSH algorithm utilizing κTDC is subsequently referred to as κTSH in literature, 30 however, we propose to call it rather κFSSH since it is the FSSH method just supplemented with κTDC. The expression in eq 12 is not commonly used and in this manuscript we will refer to it as λTDC. A FSSH algorithm utilizing λTDC will be analogously called λFSSH. As mentioned in Section 2.2, hopping algorithms based solely on energies suffer from discontinuities in PESs. While spurious minima can be detected for LZSH and hops blocked, the situation is more complicated for both κand λFSSH. κ/ λTDC requires numerical second derivatives of Zik which are strongly affected by discontinuities (usually overestimating couplings). Since TDSE is propagated in κ/λFSSH using the κ/λTDC, the discontinuities will impact the electronic coefficients. However, it is not straightforward how to mitigate such problems since the hop can happen even long after encountering the discontinuity. Although we are aware of ad hoc corrections that attempt to address this issue, 51 these solutions are typically system-specific and not generally applicable. We therefore do not include them in our simulations, consistent with ref 50, which also benchmarks the κFSSH method without any patches. As shown in the Supporting Information, the effect of the patch applied to LZSH is small and certainly does not account for the discrepancies between κFSSH and LZSH observed in azobenzene. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9787 3. COMPUTATIONAL DETAILS 3.1. Analytic Potentials. We have used several analytical potentials to showcase the performance of the investigated TSH algorithms. The simplest and also the most popular twostate system one is the original simple avoided crossing (SAC) potential developed by Tully. 16 The general diabatic potential, defined in atomic units, has the form = =V x V x A x e( ) ( ) sgn( )(1 ) x Bx 00 11 sgn( ) (14) = =V x V x Ce( ) ( ) x 01 10 2 (15) where we have used A= 0.01 and B=−1.6. The parameter C will be specified at the corresponding places, as it is used to analyze the effect of the TDC width. To further challenge TSH algorithms, we present their performance on a system featuring three-state crossing. For direct comparison with Tully’s SAC model, we have expanded it by adding a zero-energy state between Tully’s original curves (see the first column of Figure 2). The resulting potential takes the form = = =V x V x A x e V( ) ( ) sgn( )(1 ), 0 x Bx 00 22 sgn( ) 11 (16) = = = = = = V x V x V x V x Ce V V ( ) ( ) ( ) ( ) , 0 x 01 10 12 21 02 20 2 (17) where we have used A= 0.01 and B=−1.6. To ensure the corresponding three-state adiabatic potential has a zero-energy state, we set V02 =V20 = 0. Again, the coupling constant Cwill be specified at the appropriate places. We also analyzed the SH algorithms on potentials with trivial crossings. The two-state model used for the analysis has the form = + =V x x V x x( ) 0.0025( 1) , ( ) 0.0025( 1) 00 2 11 2 (18) = =V x V x( ) ( ) 0 01 10 (19) For a three-state model with trivial crossings, we used = = + = + V x x V x A x V x B x ( ) 0.0025 , ( ) 0.0025( 2) , ( ) 0.0025 00 2 11 2 22 2 (20) = = = = = =V V Ce V V Ce V V, , 0 x x 01 10 ( 2) 12 21 02 20 2 2 (21) where A,Band Care constants, which will be specified at the corresponding places. 3.2. Vibronic Coupling Model. The vibronic coupling (VC) model is an analytical Hamiltonian, typically formulated in the diabatic basis, that is fitted to describe the vibrational motion of a system. 59−61 The model allows for highly efficient numerical evaluation of electronic structure quantities, even for highly dimensional systems. In this work, we employ a VC Hamiltonian representing an 8-dimensional model of the uracil cation. In the VC approach, the Hamiltonian is generally expressed as =++++µH H W W W (0) (1) (2) (3) (22) where H(0) is the reference potential and the Wterms represent the various vibronic coupling contributions. 60 The Hamiltonian is typically written in terms of the mass-frequency scaled coordinate vector Q. We can write the reference potential as = +H E VQ Q( ) ( ) jj j l l j( ) (23) Figure 2. Nonadiabatic dynamics simulations on Tully’s Simple Avoided Crossing potential defined in (14) using two different sets of diabatic coupling constants for eq 15. The first row uses C= 0.005 a.u. while the second row uses C= 0.0005 a.u. The first column displays the potential energies in the adiabatic basis. The second column presents the TDC for a selected trajectory, while the third column shows the ground-state population throughout the simulation. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9788 with Ejbeing the energy shift of an electronic state jand Vl, depending on the anharmonicity of the mode l, being either harmonic or Morse potential written as = = [ ] +V Q V d e 1 2e 1 l j lll j l ja Q Q l j( ) 2( ) ( ) ( ) 2 ( ) l j l l j ( ) ,0 ( ) (24) where ωlis the vibrational frequency of mode land dl (j),al (j), Ql, 0 (j),el (j)are constants specific for the vibrational mode land electronic state j. The first-order coupling terms W(1) are given by =W Q jj l l j l (1) ( ) (25) =W Q jk l l jk l (1) ( ) (26) where κl (j)is a constant for the vibrational mode lin an electronic state jand λl (jk)denoting a coupling constant between electronic states jand kfor a vibrational mode l. If the Hamiltonian in eq 22 is truncated after the first-order terms, the model is referred to as the linear vibronic coupling (LVC) model. 62 In this work, we also consider diagonal elements of W(2) and W(4) in the form of =W Q 1 2 jj l l j l (2) ( ) 2 (27) =W k Q 1 24 jj l l j l (4) ( ) 2 (28) where γland μlare constants. All modes, constants, and reference potentials used in the VC models in this work are provided in the Supporting Information. 3.3. Nonadiabatic Dynamics on Analytic Potentials. We performed numerically exact nonadiabatic dynamics by solving TDSE in the diabatic basis using the split-operator method. 63,64 In all simulations, the initial wave function was chosen as a Gaussian wavepacket in a form = +x t x x ip x x( , ) exp( ( ) ( )) 0 0 2 00 (29) where x0and p0represent the initial position and momentum, respectively, iis the imaginary unit and Ωis the normalization factor. Simulations were carried out on a spatial domain of [−24 a.u., 48 a.u.] (or [−10 a.u., 10 a.u.] for trivial crossings) discretized into 8192 grid points, using a mass of 2000 a.u. and a time step of 10 a.u. For Tully’s potentials (eqs 14 and 16) we have set the initial momentum to 15 a.u., all other simulations on model potentials used zero initial momentum. To carry out FSSH simulations on analytic potentials, we employed the norm preserving interpolation (NPI) scheme to approximate the TDC. 65,66 Adiabatic surface gradients were obtained using first-order finite differences with a step size of 0.001 a.u. The time step used was set to 1 a.u. and each simulation on analytic potential used 10,000 trajectories. 3.4. Nonadiabatic Molecular Dynamics on Ab Initio Potentials. Three molecular systems are investigated with ab initio nonadiabatic dynamical simulations in this work: photoisomerization of S1excited cis-stilbene and transazobenzene, and photodissociation of S2excited cyclobutanone. For cis-stilbene. We employed the multireference configuration interaction with singles and doubles (MRCISD) method based on the semiempirical OM3 Hamiltonian and state-averaged complete active space self-consistent field (SACASSCF) approach. The SA-CASSCF calculations utilized an active space comprising two electrons in two πand π*orbitals, denoted as (2,2), reproduced from a previous study. 57 In contrast, the OM3-MRCISD method, due to its computational efficiency, enabled the use of a significantly larger active space, such as (12,17), also reported in earlier work. 39 All cis-stilbene trajectories were initiated in the S1state (considering only the lowest two states), with initial geometries and momenta sampled from a harmonic Wigner distribution around the MP2/6-31G*minimum of cis-stilbene, assuming a temperature of 298.15 K, see ref 38, for more details. The SA-CASSCF calculations were carried out using the BAGEL package, 67 while OM3-MRCISD computations were performed with the MNDO v7.0 code. 68 All dynamical simulations were executed within our molecular dynamics code ABIN and the time step was set to 5 a.u. 69 For the trans-azobenzene, we have performed only the OM3-MRCISD calculations, also with the active space of (12,17). 39 All the trajectories start in the S1state (considering only the lowest two states) and the initial conditions are sampled using a harmonic Wigner distribution around the B3LYP/6-31G*minimum of trans-azobenzene, also assuming a temperature of 298.15 K. The cyclobutanone simulations with extended multistate complete active space second-order perturbation theory (XMS-CASPT2) were done as described in ref 37. For the CASSCF simulations, we employed the same active space, initial conditions, and general strategy as in ref 37, but we did not restart the trajectories in the ground state. Each OM3-MRCISD simulation was performed using 500 trajectories. The number of failed trajectories for these simulations did not exceed 20, and all of them were due to poor energy conservation. The cis-stilbene SA-CASSCF simulations employed 168 FSSH trajectories, 89 κTSH trajectories, and 91 LZSH trajectories. All cyclobutanone simulations comprised 119 trajectories. 3.5. Trajectory Analysis. The principal quantities of interest in these simulations are the electronic state populations. At any given time t, the population of a specific state is computed by taking the fraction of trajectories in that state relative to the total number of trajectories still active. In the event that a simulation terminates immediately after a jump to the ground state, we assume that the trajectory remains in the ground state for the rest of the simulation. Consequently, this trajectory continues to contribute to the ground-state population even after the simulation has failed. Because each population can be viewed as a realization of a multinomial random variable at each time point, we estimate its confidence interval by invoking the normal approximation. Specifically, for the kth state, the confidence interval at time tis given by = ±p t z p t p t N t ( ) ( )(1 ( )) ( ) k k k (30) where pk(t) denotes the population of state kat time t,N(t) is the number of trajectories still considered at the time t, and zis the 97.5% quantile of the standard normal distribution (i.e., z ≈1.96) for a two-sided 95% confidence interval. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9789 4. RESULTS We begin by applying the LZSH, FSSH, and κFSSH algorithms to various analytic potentials, evaluating their theoretical performance against the numerically exact solution of the time-dependent Schrodinger equation. Next, we extend our analysis to more realistic systems, examining how the curvature-based algorithms perform on the LVC model of Uracil. Finally, we shift to real molecular systems which feature less forgiving PESs for the curvature-based algorithms. Finally, the comparison of these algorithms in ab initio simulations of stilbene, azobenzene, and cyclobutanone will be provided. 4.1. Model Potentials. 4.1.1. Two-State Avoided Crossing Model. Let us start by examining how the various algorithms perform on a two-level system, specifically Tully’s SAC potential (see Section 3.1). The results for this two-level case with diabatic coupling constant C= 0.005 a.u. are shown in the top row of Figure 2. In the bottom row of the same figure, we use the same potential expression but with coupling constant C= 0.0005 a.u. This allows us to assess how sensitive the final electronic populations are to the exact shape of the TDC. Inspecting the TDCs from the top row of Figure 2, it is clear that NPI TDC is wider than the κTDC. This discrepancy directly affects the ground-state populations (right panels) and κFSSH performs poorly compared to LZSH and especially FSSH. To confirm that the culprit is indeed the TDC approximation, we repeated the simulations with a modified diabatic coupling which causes the TDC to be more time localized. This improves the agreement between κTDC and NPI TDC simulations, shown also in ground-state populations although one still observes that the LZSH algorithm outperforms the κFSSH method. Note that λFSSH fully correlates with κFSSH. Hence, LZSH appears to be a more robust choice in this particular two-state problem. In both cases, FSSH with NPI TDC aligns with the exact results. 4.1.2. Three-State Avoided Crossing Model. Next, we push the algorithms further by considering a three-state analogue of Tully’s SAC potential, defined in eqs 16 and 17 with C= 0.003 a.u. The corresponding results are shown in Figure 3. As with the two-level model, we also present (in the bottom row) results obtained by using a modified coupling with C= 0.0003 a.u. Once again, κTDC is too narrow compared to NPI TDC. However, we now encounter additional problems stemming from the three-state nature of the model. While the NPI method correctly predicts TDC to be zero between states 0 and 2, the κapproach assigns nonzero TDC elements to every pair of states. A closer inspection of the underlying potential clarifies these failures. Analytically, one can show that for the chosen threestate model, the NACV between states 0 and 2 is vanishing. The NPI method captures this fact; the κTDC approach, based solely on adiabatic curvature, does not. One can further demonstrate that = = Z Z Z Z Z Z 01 01 02 02 12 12 . Under these conditions, the κTDC is the same for all pairs of states. We should also note that the κTDC underestimates the magnitude of TDC, which will become a problem in the following analysis. Finally, we add that λTDC and κTDC again exhibited the same shape. The results also highlight an excellent performance of LZSH in its maximum second derivative version (LZSHmaxcurv), contrasting with the poor performance of κFSSH. The success of LZSHmaxcurv dwells in its ability to recognize the interacting states by maximum curvature of the energy difference. The nearest state version introduced in Section 2.2 performs poorly in this case since it transfers population solely to state 1 and neglects state 2 completely. Similarly, we have tested a version where the hop goes to a state with maximum hopping probability, leading again to the population of only state 1. In both cases, FSSH with NPI TDC is again closest to the exact results. Hence, if analytic couplings are available, FSSH is Figure 3. Nonadiabatic dynamics simulations on a three-state crossing potential defined in (16) using two different sets of diabatic coupling constants for eq 17. First row uses C= 0.003 a.u. while the second row uses C= 0.0003 a.u. The first column displays potential energy curves in the adiabatic basis. The second column presents the TDC for a selected trajectory, while the third column shows the ground-state population over the course of the simulation. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9790 the best approach from those selected. However, if couplings are not available, it appears beneficial to prioritize LZSHmaxcurv over κFSSH. 4.1.3. Trivial Crossing Models. The varying TDC magnitude across different methods naturally raises the question of how surface hopping algorithms perform when applied to potential energy surfaces exhibiting trivial crossings. In such cases, the associated nonadiabatic coupling becomes singular at the crossing point. Figure 4 presents the nonadiabatic dynamics computed for the two-state model defined in eq 18 and three-state potential defined in eq 20, using three different sets of parameters. The first row of Figure 4 corresponds to the potential defined in eq 18, which exhibits a TDC singularity at the point x= 0. As we see, the magnitude of the NPI TDC is sufficient to allow FSSH a 100% transfer to the ground state. The LZSHmaxcurv algorithm also predicts a full transfer to the ground state due to the infinitely small gap between the states at the crossing. However, as we have already seen in the Figure 4. Nonadiabatic dynamics computed on the model potential of eq 18 (top row) and eq 20 (remaining rows). The left column depicts the adiabatic potential energy curves, the middle column is the TDC of a random trajectory near the first crossing, and the right column displays the corresponding time-dependent adiabatic state populations. In the second row, the diabatic potential parameters are set to A= 0.01 a.u., B= 0.02 a.u. and C= 0 a.u. The third row uses A= 0.05 a.u., B= 0.1 a.u. and C= 0 a.u., increasing the energy separation between states. The last row introduces a small Gaussian coupling with A= 0.01 a.u., B= 0.02 a.u. and C= 2 ×10−4a.u., which perturbs the crossing topology while maintaining full deexcitation of the initially populated state. All simulations started at the highest adiabatic state at the point x=−7 a.u. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9791 Figures 2 and 3,κFSSH and λFSSH underestimate the magnitude of the TDC and the size of the coupling is simply not high enough to predict the full population transfer. Due to the insufficient magnitude of the TDC in these methods, the electronic coefficients also do not reach the full transfer at the first crossing, making the subsequent jumps even further from the exact solution. We should also note that the κTDC has a small tail due to the =0 Z t d d jk approximation, contrary to the λTDC. This small tail partially compensates for the insufficient peak magnitude of the coupling and enhances the performance over λTDC. However, this is only a result of the compensation of errors. In the second row of Figure 4, the potential energy surfaces exhibit two trivial crossings. The NPI TDC and LZSHmaxcurv algorithm, again, reproduce the population transfer perfectly. Both κFSSH and λFSSH approaches perform poorly in this scenario, owing mainly to the size of the coupling, leaving residual population on S2rather than effecting a full transition. We should also note here that κTDC and λTDC predict nonzero coupling between S0and S2at the first crossing that is not high enough to allow for a population transfer to S0but affects the overall electron dynamics. In the third row of Figure 4, the same diabatic potential as in the second row is used, but the energy gap between adjacent surfaces has been increased by a factor of 5, which shifts the first crossing point away from the S2minimum. Nonetheless, since Z t d d ik 2 2 is not necessarily zero away from the crossing, the κTDC has a nonzero tail which is, in this case, not negligible. Due to this tail, we see more accurate population transfer for κFSSH than λFSSH. Still, the performance of FSSH and LZSHmaxcurv remains unchallenged. Histogram of the hopping positions around the first crossing for this model can be found in the Supporting Information. In the bottom row of Figure 4, we introduce a localized Gaussian diabatic coupling with a peak magnitude of 10−4a.u. This modification perturbs the energy crossing only slightly while preserving the full transfer at the crossing points. FSSH with NPI TDC provides, as expected, exact populations. The LZSHmaxcurv algorithm also yields quantitatively correct populations. Moreover, the Gaussian coupling renders the κTDC/λTDC profile more similar to the NPI TDC, such that κFSSH/λFSSH results become acceptable. We should still note a nonzero κTDC and λTDC between the S0and S2state at the crossings, which causes small deviations of populations using these methods. Note that for all the three-state models in Figure 4, both LZSHmaxcurv and LZSHnearest perform equally since the interacting state with the largest Zik is always the nearest state. The difference in the algorithms appears only when interaction over multiple states occurs, which is not a common case in photodynamics. In summary, Figures 2 and 3illustrate the importance of obtaining a reliable TDC. The limited spatial width inherent to κTDC does introduce errors, yet the more serious failures of the κTDC and λTDC schemes arise from structural flaws that assign nonadiabatic coupling even between electronic states that should not be coupled. Figure 4 further demonstrates that, in the case of trivial crossings, the κFSSH and λFSSH methods inherently fail to predict correct population transfer, whereas the LZSHmaxcurv algorithm remains robust and therefore constitutes a superior choice for such systems. Finally, we note that there is currently an initiative to apply κFSSH also for extended system 70 which, however, often exhibit trivial crossings. Vogt et al. compare FSSH results using κTDC (coined TDBA) as well as standard TDC in ref 70, showing notable discrepancies for κTDC. The authors suggest caution when using κTDC as its two state nature leads to population leakage to other states. In light of our results, we Figure 5. Population dynamics of the uracil cation computed using an eight-mode vibronic coupling model. Initial conditions were sampled from a Wigner distribution. Trajectories were propagated with a time step of 10 au for 250 iterations, following the protocol outlined in ref 71. The initial electronic state was a mixed state composed of 94% D2, 5% D1, and 1% D0populations. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.5c01176 J. Chem. Theory Comput. 2025, 21, 9784−9798 9792