We perform on-the-fly non-adiabatic molecular dynamics simulations using the symmetrical quasi-classical (SQC) approach with the recently suggested molecular Tully models: ethylene and fulvene. We attempt to provide benchmarks of the SQC methods using both the square and triangle windowing schemes as well as the recently proposed electronic zero-point-energy correction scheme (the so-called γ correction). We use the quasi-diabatic propagation scheme to directly interface the diabatic SQC methods with adiabatic electronic structure calculations. Our results showcase the drastic improvement of the accuracy by using the trajectory-adjusted γ-corrections, which outperform the widely used trajectory surface hopping method with decoherence corrections. These calculations provide useful and non-trivial tests to systematically investigate the numerical performance of various diabatic quantum dynamics approaches, going beyond simple diabatic model systems that have been used as the major workhorse in the quantum dynamics field. At the same time, these available benchmark studies will also likely foster the development of new quantum dynamics approaches based on these techniques.

Simulating on-the-fly non-adiabatic quantum dynamics in molecular systems remains a central challenge in modern theoretical chemistry despite the impressive progress made in the past several decades.1–28 The two main components for performing an on-the-fly quantum dynamics simulation are (i) obtaining accurate electronic structure information and (ii) using it to propagate the coupled motion of nuclear and electronic degrees of freedom (DOF) in an efficient manner.29 Mixed quantum classical (MQC) approaches, such as the fewest-switches surface hopping1 (FSSH) and the mean-field Ehrenfest30 (MFE) approaches, which use the outputs of electronic structure methods to evolve the electronic subsystem quantum mechanically and nuclear DOFs classically, have remained popular for simulating on-the-fly quantum dynamics. However, the inherent mixed-quantum classical approximation in these approaches can lead to the break-down of detail balance,31 the artificial creation of electronic coherence,18 or incorrect chemical kinetics.18 

In response to these deficiencies, a wide range of non-adiabatic dynamics approaches have been developed in the diabatic representation, some of which include the partial linearized density matrix7,32 (PLDM), symmetrical quasi-classical11,12 (SQC) approach, state-dependent ring polymer molecular dynamics,13,15,23 quantum-classical path integral (QCPI) approach,33–36 and quantum classical Liouville equation (QCLE) dynamics.14,37 In particular, the recently developed γ-SQC has been shown27 to provide impressively accurate non-adiabatic photo-dissociation quantum dynamics with coupled Morse potentials through the adjusted zero point energy parameter of the mapping variables, thus appearing to be a promising method to simulate on-the-fly quantum dynamics of complex molecular systems. Testing these approaches with simple model systems becomes the major workhorse in the quantum dynamics field. What have been largely missing, on the other hand, are the calculations that go beyond simple diabatic models. However, reformulating these approaches from diabatic to adiabatic representation often requires non-trivial theoretical tasks and introduces new numerical complications.

In our recent works, we have developed the quasi-diabatic (QD) propagation scheme38–42 as a general strategy to seamlessly combine a diabatic quantum dynamics approach with the adiabatic outputs of an electronic structure method. The QD propagation scheme uses the adiabatic states with a reference nuclear geometry (the so-called “crude adiabatic” states) as the local diabatic states during a short-time propagation and continuously updates the QD basis at each consecutive nuclear propagation step. In this propagation scheme, one does not construct a global diabatic representation but uses a sequence of local diabatic representations for each short-time segment to propagate quantum dynamics. Note that the quasi-diabatic propagation scheme38–42 should not be confused with the approximate diabatic representation, which is also often referred to as the QD representation in the literature.43–45 

In this work, we use the QD propagation scheme to combine the γ-SQC approach with adiabatic outputs of the complete active space self-consistent field (CAS-SCF) approach to perform on-the-fly non-adiabatic quantum dynamics in molecular systems. We directly simulate photo-excited non-adiabatic dynamics in two molecular systems, ethylene and fulvene. These molecular systems have been recently proposed46 to be analog of the Tully curve-crossing models,1 I and III, respectively. The original Tully models1 explore nuances in excited state dynamics that are ubiquitous in “real” molecular systems to varying degrees of complexity, while only involving one degree of freedom, and have been extensively used to benchmark quantum dynamics approaches.6,7,11,47–49 The molecular analogs of the Tully models, on the other hand, capture the basic physics of the original Tully models while concurrently showcasing complex dynamical features due to the coupled motion between the electronic and multiple nuclear DOFs. These molecular systems serve as robust benchmarks offering complex non-adiabatic dynamics beyond one-dimensional, overly simplified model systems and are representative of typical molecular systems.

Our numerical results demonstrate that the zero-point energy (ZPE) corrected SQC (γ-SQC) improves the population dynamics compared to the original SQC approach, benchmarked against the more accurate ab initio multiple spawning (AIMS) approach. Our numerical results also show that γ-SQC can outperform the state-of-the-art decoherence-corrected surface hopping (dTSH) approach, with a similar numerical cost in terms of the number of trajectories. Overall, our results demonstrate the accuracy and applicability of the γ-SQC approach for ab initio on-the-fly simulation enabled by the QD propagation scheme, opening up future opportunities for simulating on-the-fly quantum dynamics of complex molecular systems.

The SQC approach11,12 uses symmetrical window functions to sample electronic DOF at initial time and provides an estimate of the reduced density matrix at later times. It relies on the Meyer–Miller–Stock–Thoss (MMST) mapping Hamiltonian approach, which transforms the electronic degrees of freedom onto an effective set of singly excited and fictitious classical harmonic oscillators. It has been shown to provide accurate non-adiabatic dynamics in a wide range of model systems.11,50–52 Recently, the original SQC method was also combined with the QD propagation scheme41 for direct quantum dynamics simulation seamlessly using adiabatic electronic structure outputs.42 Here, we briefly discuss the essential idea of the mapping Hamiltonian, the γ-SQC approach, and the QD propagation scheme.

The Meyer–Miller–Stock–Thoss (MMST) mapping representation53–55 transforms the discrete electronic DOFs onto an effective set of fictitious, singly excited classical harmonic oscillators, thus mapping the electronic non-adiabatic dynamics onto these oscillators’ phase space motion.

The total molecular Hamiltonian in the diabatic representation is expressed as follows:

Ĥ=T̂+ijVij(R̂)|ij|,
(1)

where Vij(R̂)=i|V̂(r̂,R̂)|j are the matrix elements of the electronic Hamiltonian in the diabatic basis {|i⟩}. Using the Meyer–Miller–Stock–Thoss54–56 mapping representation, the discrete electronic states are transformed into continuous phase-space variables,

|ij|âiâj,
(2)

where âi=(q̂iip̂i)/2 and âj=(q̂j+ip̂j)/2. With this transformation, the molecular Hamiltonian in Eq. (1) is transformed into the following MMST mapping Hamiltonian:

Ĥm=T̂+12ijVij(R̂)p̂ip̂j+q̂iq̂j2γδij,
(3)

where γ = 0.5 is the ZPE for the mapping harmonic oscillators. Historically, it is recognized as the Langer correction by Meyer and Miller56 for the quasi-classical description. Note that, until Eq. (3), no approximations have been made.

In the SQC approach, instead of evolving all DOF quantum mechanically, the coupled electronic–nuclear dynamics are propagated using the following classical Hamiltonian:52 

Hm=P22M+12ijVij(R)pipj+qiqj2γδij,
(4)

where γ is viewed as a parameter52 that specifies the ZPE of the mapping oscillators.57 

Classical trajectories are generated based on Hamilton’s equations of motion,

q̇j=Hm/pj,ṗi=Hm/qi,
(5)
Ṙ=Hm/P,Ṗ=Hm/R=F,
(6)

with the nuclear force expressed as

F=12ijVij(R)pipj+qiqj2γδij.
(7)

Overall, the MMST mapping Hamiltonian provides a consistent classical footing for both electronic and nuclear DOFs, and the non-adiabatic transitions between electronic states are captured through the classical motion of the fictitious harmonic oscillators.

To sample the electronic initial condition and estimate the population, the SQC approach uses the action-angle variables, {ej, θj}, which are related to the canonical mapping variables {pj, qj} through

ej=12pj2+qj2;θj=tan1piqi,
(8)

and the inverse relations are

qj=2ejcos(θj);pj=2ejsin(θj),
(9)

where ej is a positive-definite action variable introduced by Cotton and Miller that is directly proportional to the mapping variables’ radius in action-space,27 which allows for conceptual simplification (compared to the nj = ejγ action variable used in previous works11,41,58), as it is independent of the ZPE parameter γ, which will be allowed to be state-dependent in all subsequent sections of this work.

The SQC approach allows for the population of electronic state |j⟩ to be evaluated as52 

ρjj(t)=TrRρ̂R|ii|eiĤt/|jj|eiĤt/12πN+MdτρW(P,R)Wi(e(0))Wj(e(t)),
(10)

where ρ̂(0)=|ii|ρ̂R is the initial density operator, ρW(P, R) is the Wigner density of ρ̂R operator that contains N nuclear DOFs, e={e1,e2,,eF} is the positive-definite action variable vector for F electronic states, Wi(e) = δ(ei − (1 + γ))ijδ(ejγ) are the Wigner transformed action variables,59 and dτdP · dR · de · dθ.

For practical reasons, the above delta functions are artificially broadened using two well-explored distribution functions (i.e., square and triangle) that can be used to bin the resulting electronic action variables in action-space, depicted for any two-state projection in Figs. 1(a) and 1(b).52 The square distribution for an F-state system is defined as27 

Wj(e)=w1(ej)jjFw0(ej),
(11)

where the function wN(e) is expressed as

wN(e)=1,0<eN<2×0.366,0else,
(12)

and γ=3/210.366 is the optimal width parameter of the square window.11,28

FIG. 1.

(a) Symmetric quasi-classical (a) square and (b) triangle window distributions depicted for a two-state projection involving states 1 (blue) and 2 (red). The action-space window distributions are depicted using the positive-definite action variables {ek}, which are shifted quantities according to the corresponding zero-point energy parameter γ.

FIG. 1.

(a) Symmetric quasi-classical (a) square and (b) triangle window distributions depicted for a two-state projection involving states 1 (blue) and 2 (red). The action-space window distributions are depicted using the positive-definite action variables {ek}, which are shifted quantities according to the corresponding zero-point energy parameter γ.

Close modal

The triangle window27,58 is expressed as

Wj(e)=w1(ej)jjFw0(ej,ej),
(13)

where

w1(e)=(2e)2F,1<e<2,0else
(14)

and

w0(e,e)=1,e<2e,0else,
(15)

and trajectories are assigned to state j at time t if ej ≥ 1 and ej < 1 for all j′ ≠ j. The ZPE parameter in the triangle scheme is γ = 1/3. The triangle window scheme for a two-state system is depicted in Fig. 1(b).

The time-dependent population at time t is then calculated by applying the window function estimator to action variables {ej(t)} for an ensemble of trajectories. Starting from the initial diabatic state |i⟩, the time-dependent population of the states |j⟩ is computed using Eq. (10). However, by using the window function estimator, the total population is no longer properly normalized due to the fraction of trajectories that are outside of any window region at any given time.11 Thus, the total population must be normalized11 with the following procedure:

ρjj(t)/i=1Nρii(t)ρjj(t).
(16)

This SCQ approach provides a dramatic improvement to Ehrenfest dynamics, even though they utilize the same equation of motion for the coupled electronic–nuclear DOFs.52,60 The SQC method allows for the elimination of known issues present in Ehrenfest dynamics, including preserving detailed balance,31,60,61 and this method has shown to be quite accurate in an assortment of model systems50,51,58,62 while only needing a few thousand trajectories for convergence.11,51,52 As such, the SQC method is well-suited for use in on-the-fly non-adiabatic simulations of real molecular systems.

It was recently proposed that the mapping zero-point energy should be chosen in such a way as to constrain the initial force to be composed purely from the initially occupied state,27 which was not previously enforced using a fixed γ in Eq. (7). This new scheme has shown to provide a significant improvement for photo-dissociation problems with coupled Morse potentials27 and has been combined with the kinematic momentum approach63 to carry out on-the-fly simulations of the methaniminium cation.64 

The basic logic of this scheme is to choose an appropriate γj for each state |j⟩ in a given individual trajectory, such that the initial population is forced to respect the initial electronic excitation focused onto a single excited state. If the initial electronic state is |i⟩, then

γj=ejδji,
(17)

or equivalently,

δji=ejγj,
(18)

where the {ej} are uniformly sampled, and then the γj are chosen to satisfy Eq. (18).

These γj will be explicitly used in the equations of motion (EOMs) in Eqs. (5)(7), and in particular, the nuclear forces are now

F=12ijVij(R)pipj+qiqj2γjδji,
(19)

ensuring the initial forces (at t = 0) are simply F = −∇Vii(R). Previously, without any adjustments to γk, the chosen values for γk were only dependent on the windowing function itself, i.e., γk = 0.366 for the square windows and γk = 1/3 for the triangle windows. With the above γ-correction method,27 each individual trajectory will have its own state-specific γj for state |j⟩ that is completely independent of the choice of the window function. The choices of the γ parameter for different SQC approaches are summarized in Table I.

TABLE I.

Choice of γ for state |j⟩ with initial excitation on state |i⟩ of each SQC method employed in this work.

Methodγ-correctionγj
□-SQC False 0.366 
Δ-SQC False 1/3 
□-γ-SQC True ejδji 
Δ-γ-SQC True ejδji 
Methodγ-correctionγj
□-SQC False 0.366 
Δ-SQC False 1/3 
□-γ-SQC True ejδji 
Δ-γ-SQC True ejδji 

Note that reformulating γ-SQC in the adiabatic representation [such as the kinematically transformed (KT) SQC63] has been done recently to perform on-the-fly simulations.64 However, formulating a quantum dynamics approach in the adiabatic representation introduces additional numerical issues as the molecular Hamiltonian in the adiabatic representation involves first and second derivative coupling elements that are typically sharply peaked around avoided crossings and become singular at a conical intersection (CI). In our previous work,41 we showed that the kinematically transformed SQC (KT-SQC), formulated in the adiabatic representation, may require substantially small time step to converge dynamics due to the presence of sharply peaked first derivative coupling elements, even though the second derivative coupling elements are removed through mathematical transformation in this formalism.63 To this end, we use the quasi-diabatic propagation scheme to directly interface the diabatic γ-SQC approach with adiabatic electronic structure information.

In this work, we combine the γ-SQC approach formulated in the diabatic representation with adiabatic outputs of an electronic structure approach using the QD propagation scheme.38,41,42 This scheme has been discussed in our previous work.42 Here, for the sake of completeness of the current work, we briefly summarize the basic idea of the QD propagation scheme.

Despite recent theoretical progress,65–69 strict diabatic states {|i⟩, |j⟩} are neither uniquely defined nor routinely available for “real” molecular systems. In contrast, it is convenient to obtain adiabatic states by solving the following eigenequation:

V̂(r̂,R)|Φα(R)=Eα(R)|Φα(R),
(20)

where V̂(r̂;R) is the electronic part of the molecular Hamiltonian at a nuclear configuration R and |Φα(R)⟩ is the adiabatic state, which is the eigenstate of V̂(r̂;R), with the corresponding eigenvalue Eα(R) referred to as the adiabatic potential energy.

The central idea behind the QD scheme is realizing that, to propagate quantum dynamics with diabatic dynamics approaches, one only needs locally well-defined diabatic states, and these local diabatic states can simply be adiabatic states with a reference geometry (which are commonly referred to as the crude adiabatic states). The QD scheme38 directly uses the adiabatic states associated with a reference geometry as the local diabatic states during a short-time quantum propagation and dynamically updates the definition of the QD states along the time-dependent nuclear trajectory. Between two consecutive nuclear geometries in the time-evolution, R(t0) and R(t1), one can use the adiabatic states at R(t0) as a quasi-diabatic basis to propagate the electronic subsystem to the next nuclear geometry R(t1). In particular, the basis

|Φα(R0)|Φα(R(t0))fort[t0,t1]
(21)

is taken and used to represent the potential energy V̂(r̂,R) at all electronic time steps between R(t0) and R(t1),

Vαβ(R(t))=Φα(R0)|V̂(R(t))|Φβ(R0).
(22)

At each electronic time step, an interpolation [between nuclear geometries R(t0) and R(t1)] of the potential energy operator is obtained in the basis of R(t0), implying that all interpolated steps produce an off-diagonal electronic matrix in this basis. Most importantly, note that in this quasi-diabatic basis, the derivative couplings vanish by definition.

The linear interpolation70 between Vαβ(R(t0)) and Vαβ(R(t1)) is computed as

Vαβ(R(t))=Vαβ(R0)+(tt0)(t1t0)Vαβ(R(t1))Vαβ(R0),
(23)

where Vαβ(R0)=Φα(R0)|V̂(R(t0))|Φβ(R0)=Eα(R(t0))δαβ. The matrix elements Vαβ(R(t1)) are computed as follows:

Vαβ(R(t1))=λνSαλVλν(R(t1))Sβν,
(24)

where

Vλν(R(t1))=Φλ(R(t1))|V̂(R(t1))|Φν(R(t1))=Eλ(R(t1))δλν,

and the overlap matrix between two adiabatic electronic states (at two different nuclear geometries) are Sαλ = ⟨Φα(R0)|Φλ(R(t1))⟩ and Sβν=Φν(R(t1))|Φβ(R0). These overlap matrices are computed based on the approach outlined in Ref. 71. The nuclear gradients ∇Vαβ(R(t1)) ≡ ∂Vαβ(R(t1))/R are evaluated as

Vαβ(R(t1))=Φα(R0)|V̂(R(t1))|Φβ(R0)=Φα(R0)|V̂(R(t1))|Φβ(R0)=λνSαλΦλ(R(t1))|V̂(R(t1))|Φν(R(t1))Sβν.
(25)

Since {|Φα(R0)⟩} is a diabatic basis, the nuclear gradient will bypass these terms, and the resolution of identity ∑λλ(R(t1))⟩⟨Φλ(R(t1))| = 1 was inserted, which assumes the completeness of basis at the R(t1) geometry. Note that Eq. (25) includes all nuclear gradient terms without approximation.

Once the electronic DOFs are propagated to time t1, we adopt a new basis composed of the adiabatic states at the R(t1) geometry, {|Φα(R(t1))⟩}. At this time, all quantities are transformed from {|Φα(R0)⟩} to {|Φμ(R(t1))⟩} basis using the relation

|Φμ(R(t1))=αΦα(R(t0))|Φμ(R(t1))|Φα(R(t0)).
(26)

Since the mapping relation between the physical state and the singly excited oscillator state is |Φμ(R(t1))=aμ|0=12(q̂μ+ip̂μ)|0 and |Φα(R(t0))=aα|0=12(q̂α+ip̂α)|0, the relations for the mapping variables associated with two bases are

12(q̂μ+ip̂μ)|0=αΦα(R(t0))|Φμ(R(t1))12(q̂α+ip̂α)|0.
(27)

For molecular systems, one can always find a suitable choice for the basis set in order to make ⟨Φα(R(t0))|Φμ(R(t1))⟩ real, which guarantees that the mapping variables are transformed with the same relations as the bases. Based on this relation in Eq. (27), we transform the time-dependent mapping variables between the two consecutive QD bases as follows:

αqαΦα(R(t0))|Φμ(R(t1))qμ,
(28a)
αpαΦα(R(t0))|Φμ(R(t1))pμ.
(28b)

When performing the transformation in Eqs. (26) and (28), the eigenvectors maintain their mutual orthogonality subject to a very small error when they are expressed in terms of the previous basis due to the incompleteness of the basis.8,72 Nevertheless, the orthogonality remains to be well satisfied among {|Φα(R(t0))⟩} or {|Φλ(R(t1))⟩}. This small numerical error generated from each step can, however, accumulate over many steps and cause a significant error at longer times, leading to non-unitary dynamics.8,72 This potential issue can be easily resolved by using the orthonormalization procedure among the vectors of the overlap matrix ⟨Φα(R(t0))|Φμ(R(t1))⟩, as has been done in our previous work39 for simulating photoinduced charge transfer dynamics. Here, we perform the Löwdin orthogonalization procedure73 as commonly used in the local diabatization approach72 to ensure this.

As the nuclear geometry closely follows the reference geometry throughout the propagation, the QD basis forms a convenient and compact basis. Note that, in principle, one needs infinite crude adiabatic states {|Φα(R0)⟩} to represent the time-dependent electronic wavefunction, because the electronic wavefunction could change rapidly with the motion of the nuclei, and the crude adiabatic basis is only convenient when the reference geometry R0 is close to the nuclear geometry R. By dynamically updating the basis in the QD scheme, the time-dependent electronic wavefunction is expanded with the “moving crude adiabatic basis”74 that explores the most relevant and important parts of the Hilbert space, thus requiring few states for quantum dynamics propagation.

Thus, the QD representation provides several unique advantages over the strict diabatic or adiabatic representation for quantum dynamics propagation. On one hand, the QD basis is constructed from the crude adiabatic basis, which can be easily obtained from any commonly used electronic structure calculation. On the other hand, the diabatic nature of the QD basis makes derivative couplings explicitly vanish and allows using any diabatic dynamics approaches to perform on-the-fly propagation. Furthermore, the QD scheme ensures a stable propagation of the quantum dynamics compared to directly solving it in the adiabatic representation. This is due to the fact that directly solving electronic dynamics in the adiabatic state requires the non-adiabatic coupling Φβ(R(t))|tΦα(R(t))=dβα(R)Ṙ, which might exhibit highly peaked values and cause large numerical errors20,75 when using a linear interpolation scheme.76 The QD scheme explicitly alleviates this difficulty by using the well behaved transformation matrix elements ⟨Φβ(R(t1))|Φα(R(t2))⟩ instead of Φβ(R(t))|tΦα(R(t)).

Note that the SQC approach has been derived in the adiabatic representation,63 which contains the derivative of the derivative coupling ∇dαβ(R) in the equation of motion. With kinematic momentum transformation, it explicitly eliminates the presence of the ∇dαβ(R) term in the nuclear force (instead of ignoring it), and the nuclear EOM of the KM-SQC approach is equivalent to the nuclear forces in QD-SQC.41 On the other hand, the KM-SQC EOMs still explicitly contain dβα(R) through the presence of Φβ(R(t))|tΦα(R(t)), which could lead to numerical instabilities when these derivative couplings are highly peaked. This has been extensively investigated in Ref. 41 using Tully’s avoid crossing model, which has a very narrow derivative coupling, such that it can drastically change on a time-scale that is shorter than the nuclear time step dt. When using a large dt in the KM-SQC approach, the nuclear position can encounter drastically different values of the derivative coupling from one step to another that allow for a discontinuous spike at a CI or even completely step over it,17 resulting in different long-time populations and oscillatory behavior of errors. Thus, the approaches that explicitly require derivative couplings (and those that use a simple linear interpolation scheme for obtaining them, as we have previously implemented for the KM scheme) either encounter numerical challenges or start to accumulate numerical errors.17 The QD scheme, on the other hand, provides more accurate results even when using a relatively larger time step dt, simply because the QD scheme only requires the well-behaved transformation matrix elements ⟨Φ1(R(t0))|Φ2(R(t1))⟩ instead of the highly peaked derivative coupling d12(R). That being said, there may be good alternative approaches to achieve the same attractive features for dynamics propagation, such as those recently developed norm-preserving interpolation schemes.17,20 The QD scheme is, perhaps, still the most straightforward one that allows robust dynamical propagation and enables a seamless interface between the diabatic quantum dynamics approach (such as SQC) and adiabatic electronic structure calculations.

Over the last decade, the QD representation has been used in many contexts. These schemes have been used primarily in the surface hopping approach in order to propagate the electronic coefficients.8,70,72,77,78 Furthermore, it has been also used in Gaussian wavepacket dynamics10,74,79–82 approaches where it was called the “moving crude adiabatic”74 representation. These approaches all utilize a locally diabatic basis to overcome issues arising from non-adiabatic coupling, but our approach concentrates on resolving the incompatibility between diabatic quantum dynamics approaches and adiabatic electronic structure calculations.38 

Non-adiabatic molecular dynamics simulations based on the QD-γ-SQC approach are performed using an in-house modified version42 of the Surface Hopping including Arbitrary Couplings (SHARC) non-adiabatic molecular dynamics code, interfaced to the MOLPRO electronic structure package.83,84 On-the-fly electronic structure calculations are performed at the level of the complete active space self-consistent field (CASSCF) approach, with 3SA-CASSCF(2,2)/6-31G* and 2SA-CASSCF(6,6)/6-31G* levels of theory for ethylene and fulvene, respectively.46 The CAS self-consistent calculation is performed over three lowest adiabatic states for ethylene and over two lowest adiabatic states for fulvene, whereas the quantum dynamics for both molecules are confined in the {S0, S1} subspace. All of the energies and gradients are computed at this level of electronic structure theory. The default accuracy for both the nuclear gradients and non-adiabatic vectors is 10−7 a.u.; when this criterion is not satisfied, a maximum of 900 additional wavefunction optimization iterations are used to make sure the convergence of 10−4 a.u. is reached. All of the electronic structure calculations performed during our quantum dynamics simulations converge successfully under the above criteria.

The initial Wigner distribution is sampled from the ground vibrational state ν = 0 on the ground electronic state |S0⟩, where the normal mode frequencies (in the harmonic approximation) are calculated based on the approach outlined in Refs. 85 and 86, as implemented in the SHARC package. The normal mode frequencies are computed at the level of MP2/6-31++G** with the MOLPRO package, with the optimized structure obtained at the same level of electronic structure theory for the ground state |S0(R)⟩. In particular, the nuclear density ρW(R̃,P̃) in terms of the molecular normal-mode frequencies {ω̃k} and phase space variables {R̃,P̃} is given as87 

ρW(R̃,P̃)k=1Nexptanhβω̃k2mω̃kR̃k2+1mkω̃kP̃k2.
(29)

The initial distribution {R, P} is then obtained by transforming {R̃,P̃} from the normal mode representation to the primitive coordinates using the unitary transformation that diagonalizes the Hessian matrix. The coordinates and momenta for ethylene were sampled from the above Wigner distribution of these normal modes. For fulvene, only the coordinates were sampled and the momenta were set to zero, such that the molecule encounters only a slanted conical intersection.46 

A total number of 500 trajectories are used in the QD-γ-SQC simulations to achieve converged population dynamics. A rough convergence of the population dynamics can be already achieved within 50–100 trajectories for both molecules. The nuclear time step used in the QD-γ-SQC is dt = 0.1 fs, with 200 electronic time steps for the mapping variables’ integration during each nuclear time step. The electronic structure calculations are performed only at the nuclear time step. Additional analysis, including the fraction of the trajectories outside the windows, as well as the energy conservation analysis, is provided in the supplementary material.

The overlap matrix of CAS wavefunctions between two successive nuclear time steps is calculated by using the approach outlined in Ref. 71, as implemented in the program wfoverlap. The random phases generated from electronic structure calculations for the eigenfunctions are carefully calculated and accumulated. To ensure the orthonormalization among the vectors of the overlap matrix ⟨Φα(R(t0))|Φμ(R(t1))⟩, we perform the Löwdin orthogonalization.72 All of the above routines were used as implemented in the SHARC program.88,89

In this work, we chose two molecular systems, (i) ethylene and (ii) fulvene, that were previously investigated as the ab initio analogies to Tully’s curve crossing models,46 as illustrated in Figs. 2(a) and 2(b). In panel (c), the time-dependent adiabatic potential energies of the ethylene molecule along a single nuclear trajectory for S0 and S1 states are presented. In the top sub-panel, the excited state population (blue) computed with the QD-γ-SQC approach is plotted, which shows an oscillation of the population in the avoided crossing region that eventually relaxes down to the ground state. In panel (d), the time-dependent adiabatic potential energies of the fulvene molecule are shown and they exhibit many instances of potential energy barriers within regions of near degeneracy between states, leading to reflection given a sufficiently low momentum, similar to the Tully’s model III (which is described canonically by a single spatial coordinate).1 The previously encountered avoided crossing is then visited again where the population transfers back into the excited state. This occurs twice (within the allotted time-frame) in the single trajectory presented in panel (d), which gives rise to the important features found in the population dynamics discussed in more detail later. These two molecules provide rich physics of two vastly different levels of dynamical complexity. It is important to note that the features that the Tully models present are with respect to a single nuclear coordinate, but in the ab initio molecular models, the analogy is represented through the time-dependent adiabatic potential energies, demonstrated in Figs. 2(c) and 2(d).

FIG. 2.

Adiabatic potentials for (a) Tully’s model I (a single avoided crossing) and (b) Tully’s model III (an extended region of coupling with a reflective barrier). The molecular structures of the ab initio Tully models are depicted in the insets. Along a single QD-γ-SQC trajectory, the population of the S1 (blue) state and the energies of the S0 (black) and S1 (red) states as functions of time for (c) ethylene and (d) fulvene are presented.

FIG. 2.

Adiabatic potentials for (a) Tully’s model I (a single avoided crossing) and (b) Tully’s model III (an extended region of coupling with a reflective barrier). The molecular structures of the ab initio Tully models are depicted in the insets. Along a single QD-γ-SQC trajectory, the population of the S1 (blue) state and the energies of the S0 (black) and S1 (red) states as functions of time for (c) ethylene and (d) fulvene are presented.

Close modal

Figure 3 presents the adiabatic population dynamics of ethylene upon photo-excitation to the S1 state, obtained from various SQC approaches (black solid lines), compared to the ab initio multiple spawning (AIMS) (thick green lines), an approximate Gaussian wavepacket-based non-adiabatic method, which is used as the benchmark result. The decoherence-corrected trajectory surface hopping (dTSH) approach (blue lines) is also presented for comparison, which is obtained with 660 trajectories.46 Both AIMS and dTSH results were directly adapted from Ref. 46.

FIG. 3.

Population dynamics of the S1 state in ethylene, using the square [(a) and (b)] and triangle [(c) and (d)] windowing schemes. Panels (a) and (c) utilize the original SQC method with fixed zero-point energy (ZPE) parameter γ, while (b) and (c) are computed using the trajectory-adjusted ZPE parameter γ. Note that the AIMS and dTSH data are adapted directly from Ref. 46.

FIG. 3.

Population dynamics of the S1 state in ethylene, using the square [(a) and (b)] and triangle [(c) and (d)] windowing schemes. Panels (a) and (c) utilize the original SQC method with fixed zero-point energy (ZPE) parameter γ, while (b) and (c) are computed using the trajectory-adjusted ZPE parameter γ. Note that the AIMS and dTSH data are adapted directly from Ref. 46.

Close modal

Figure 3(a) presents the results obtained from the original11,12 SQC approach using the square window scheme, where the mapping ZPE γ = 0.366 is kept as a constant for all states and trajectories. Figure 3(b) presents the results obtained from the SQC approach using the triangle window scheme58 with a fixed γ = 1/3. Both approaches provide reasonably accurate non-adiabatic dynamics compared to the AIMS results, as well as to the dTSH simulations. In particular, the □-SQC method seems to show an increased relaxation time compared to AIMS, whereas the triangle window scheme presented in panel (b) is more accurate for simulating the ethylene non-adiabatic dynamics, compared to the square window scheme presented in panel (a). This trend is in agreement with the empirical results of the recent numerical tests of both window schemes with a wide range of diabatic models,58,90 especially for models with weak non-adiabatic coupling.58 

Figures 3(c) and 3(d) present the ZPE-corrected QD-γ-SQC dynamics, obtained with the square window [panel (c)] and the triangle window [panel (d)]. For the square window scheme, we find that the γ-SQC approach [panel (c)] provides much better agreement with the AIMS benchmark compared to the original SQC method [panel (a)]. For the triangle window scheme, we find very similar short-time relaxation curves for the Δ-SQC [panel (b)] method compared to Δ-γ-SQC [panel (d)]. For the particular case of the ethylene photo-dissociation dynamics, the ZPE correction does not further improve the results when the triangle window is used, similar to the recent work that utilized the kinematic momentum formulation of SQC.64 With both Δ-SQC [panel (b)] and Δ-γ-SQC [panel (d)], we see a near quantitative agreement with AIMS, even slightly out-performing the commonly used dTSH (blue), which is obtained with 180 trajectories.46 

Figure 4 presents the non-adiabatic photo-relaxation dynamics of fulvene, which has recently been proposed as a molecular example of Tully’s model III.46 In fulvene, there exists a so-called slanted CI, whereby the wavepacket becomes reflected and re-interacts with this CI at nearly periodic times later (∼20-fs intervals), leading to the breakdown of the mixed quantum-classical methodologies due to wavepacket bifurcation as well as the added effects of encircling the CI. These non-adiabatic methods assume a single Gaussian wavepacket basis for describing the nuclear motion – higher-order modes or multiple Gaussian wavepackets are needed to fundamentally capture these effects. For the single SQC trajectory presented in Fig. 2(d), one can see that the population (blue solid line) resides in the ground state while traversing the CI but jumps back to the excited state after the interaction. Eventually, the CI will lead to the permanent relaxation to the ground state.

FIG. 4.

S1 Population dynamics for fulvene using the square [(a) and (b)] and triangle [(c) and (d)] windowing schemes. Panels (a) and (c) utilize the original SQC method with fixed zero-point energy (ZPE) parameter γ, while (b) and (c) are computed using the trajectory-adjusted ZPE parameter γ. Note that the AIMS and dTSH data are adapted directly from Ref. 46.

FIG. 4.

S1 Population dynamics for fulvene using the square [(a) and (b)] and triangle [(c) and (d)] windowing schemes. Panels (a) and (c) utilize the original SQC method with fixed zero-point energy (ZPE) parameter γ, while (b) and (c) are computed using the trajectory-adjusted ZPE parameter γ. Note that the AIMS and dTSH data are adapted directly from Ref. 46.

Close modal

Figure 4(a) presents the results obtained from □-SQC, and Fig. 4(b) presents the results from Δ-SQC. Both are providing accurate dynamics at the short time compared to AIMS, including the small shoulder in S1 population at t ∼ 10 fs as well as the subsequent plateaus in S1 population. dTSH results (blue) also provide a reasonable description for the short-time dynamics. However, at a longer time for t > 10 fs, both SQC and dTSH deviate from the AIMS results, although the SQC approach (both the square and the triangle windows) outperform dTSH. Even using different decoherence schemes in dTSH seems unable to provide further improvement, as shown in Ref. 91.

Figures 4(c) and 4(d) present the SQC dynamics with the γ correction. Contrary to the ethylene case, the SQC dynamics is significantly improved upon the γ-correction for both window schemes. In particular, the □-γ-SQC method [panel (c)] provides the most quantitative accuracy among all of the SQC-related methods, capturing both the initial inversion of population as well as the population plateau around 20 fs and finally the tail plateau around 40 fs, slightly outperforming the Δ-γ-SQC method [panel (d)]. Both γ-SQC methods greatly outperform the dTSH simulation, as there is a severe underestimation of the population transfer in the dTSH method.46 

In this work, we use the quasi-diabatic propagation scheme to directly interface the diabatic symmetric quasi-classical (SQC) approach with the electronic zero-point energy correction (the γ correction)27 and the CASSCF on-the-fly electronic structure calculations to propagate ab initio non-adiabatic dynamics. We have performed simulations for two recently suggested molecular models, ethylene and fulvene, that are closely related to the well-known simple curve crossing models of Tully. We have shown that the γ-SQC method based on the trajectory-adjusted electronic zero-point energy in classical Meyer–Miller vibronic dynamics provides very accurate non-adiabatic population dynamics when compared to ab initio multiple spawning (AIMS) and even outperforms the widely used adiabatic decoherence-adjusted surface hopping (dTSH) method. Specifically, for the fulvene molecule (which is a molecular analog of Tully’s model III), we found that the γ-correction significantly improved the accuracy of the original SQC approach for both the square and triangle window schemes. These calculations provide useful and non-trivial tests to systematically investigate the numerical performance of various diabatic quantum dynamics approaches, going beyond the simple diabatic model systems that have historically been the major workhorse in the quantum dynamics field. At the same time, these available benchmark studies will also likely foster the development of new quantum dynamics approaches.

See the supplementary material for additional analysis of the SQC trajectories and analysis of the energy conservation.

This work was supported by the National Science Foundation CAREER Award under Grant No. CHE-1845747 as well as by a Cottrell Scholar Award (a program by Research Corporation for Science Advancement). B.M.W. acknowledges support from the graduate fellowship of the Department of Physics and Astronomy at the University of Rochester. A.M. acknowledges support from the ACS Chemical Computing Group Excellence Award for Graduate Students. P.H. acknowledges support from the ACS COMP OpenEye Outstanding Junior Faculty Award. Computing resources were provided by the Center for Integrated Research Computing (CIRC) at the University of Rochester.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
2.
M.
Ben-Nun
,
J.
Quenneville
, and
T. J.
Martínez
,
J. Phys. Chem. A
104
,
5161
(
2000
).
3.
D. A.
Micha
,
J. Phys. Chem. A
103
,
7562
(
1999
).
4.
S.
Bonella
and
D. F.
Coker
,
J. Chem. Phys.
114
,
7778
(
2001
).
5.
G. A.
Worth
,
M. A.
Robb
, and
I.
Burghardt
,
Faraday Discuss.
127
,
307
(
2004
).
6.
N.
Ananth
,
C.
Venkataraman
, and
W. H.
Miller
,
J. Chem. Phys.
127
,
084114
(
2007
).
7.
P.
Huo
and
D. F.
Coker
,
J. Chem. Phys.
135
,
201101
(
2011
).
8.
F.
Plasser
,
G.
Granucci
,
J.
Pittner
,
M.
Barbatti
,
M.
Persico
, and
H.
Lischka
,
J. Chem. Phys.
137
,
22A514
(
2012
).
9.
K.
Saita
and
D. V.
Shalashilin
,
J. Chem. Phys.
137
,
22A506
(
2012
).
10.
D. V.
Makhov
,
W. J.
Glover
,
T. J.
Martinez
, and
D. V.
Shalashilin
,
J. Chem. Phys.
141
,
054110
(
2014
).
11.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
139
,
234112
(
2013
).
12.
S. J.
Cotton
and
W. H.
Miller
,
J. Phys. Chem. A
117
,
7190
(
2013
).
13.
N.
Ananth
,
J. Chem. Phys.
139
,
124102
(
2013
).
14.
C.-Y.
Hsieh
and
R.
Kapral
,
J. Chem. Phys.
138
,
134110
(
2013
).
15.
J. O.
Richardson
and
M.
Thoss
,
J. Chem. Phys.
139
,
031102
(
2013
).
16.
A. R.
Menzeleev
,
F.
Bell
, and
T. F.
Miller
,
J. Chem. Phys.
140
,
064103
(
2014
).
17.
G. A.
Meek
and
B. G.
Levine
,
J. Phys. Chem. Lett.
5
,
2351
(
2014
).
18.
J. E.
Subotnik
,
A.
Jain
,
B.
Landry
,
A.
Petit
,
W.
Ouyang
, and
N.
Bellonzi
,
Annu. Rev. Phys. Chem.
67
,
387
(
2016
).
19.
T.
Nelson
,
A.
Naumov
,
S.
Fernandez-Alberti
, and
S.
Tretiak
,
Chem. Phys.
481
,
84
(
2016
).
20.
A.
Jain
,
E.
Alguire
, and
J. E.
Subotnik
,
J. Chem. Theory Comput.
12
,
5256
(
2016
).
21.
P. L.
Walters
and
N.
Makri
,
J. Chem. Phys.
144
,
044108
(
2016
).
22.
S.
Pal
,
D. J.
Trivedi
,
A. V.
Akimov
,
B.
Aradi
,
T.
Frauenheim
, and
O. V.
Prezhdo
,
J. Chem. Theory Comput.
12
,
1436
(
2016
).
23.
S. N.
Chowdhury
and
P.
Huo
,
J. Chem. Phys.
147
,
214109
(
2017
).
24.
B. F. E.
Curchod
and
T. J.
Martínez
,
Chem. Rev.
118
,
3305
(
2018
).
25.
R.
Crespo-Otero
and
M.
Barbatti
,
Chem. Rev.
118
,
7026
(
2018
).
26.
E.
Mulvihill
,
A.
Schubert
,
X.
Sun
,
B. D.
Dunietz
, and
E.
Geva
,
J. Chem. Phys.
150
,
034101
(
2019
).
27.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
150
,
194110
(
2019
).
28.
J. E.
Runeson
and
J. O.
Richardson
,
J. Chem. Phys.
152
,
084110
(
2020
).
29.
J. C.
Tully
,
J. Chem. Phys.
137
,
22A301
(
2012
).
30.
31.
P. V.
Parandekar
and
J. C.
Tully
,
J. Chem. Theory Comput.
2
,
229
(
2006
).
32.
P.
Huo
and
D. F.
Coker
,
Mol. Phys.
110
,
1035
(
2012
).
33.
R.
Lambert
and
N.
Makri
,
J. Chem. Phys.
137
,
22A552
(
2012
).
34.
R.
Lambert
and
N.
Makri
,
J. Chem. Phys.
137
,
22A553
(
2012
).
35.
T.
Banerjee
and
N.
Makri
,
J. Phys. Chem. B
117
,
13357
(
2013
).
36.
N.
Makri
,
Int. J. Quantum Chem.
115
,
1209
(
2015
).
37.
C.-Y.
Hsieh
and
R.
Kapral
,
J. Chem. Phys.
137
,
22A507
(
2012
).
38.
A.
Mandal
,
S. S.
Yamijala
, and
P.
Huo
,
J. Chem. Theory Comput.
14
,
1828
(
2018
).
39.
A.
Mandal
,
F. A.
Shakib
, and
P.
Huo
,
J. Chem. Phys.
148
,
244102
(
2018
).
40.
A.
Mandal
,
J. S.
Sandoval C
,
F. A.
Shakib
, and
P.
Huo
,
J. Phys. Chem. A
123
,
2470
(
2019
).
41.
J. S.
Sandoval C
,
A.
Mandal
, and
P.
Huo
,
J. Chem. Phys.
149
,
044115
(
2018
).
42.
W.
Zhou
,
A.
Mandal
, and
P.
Huo
,
J. Phys. Chem. Lett.
10
,
7062
(
2019
).
43.
Y.
Guan
,
H.
Guo
, and
D. R.
Yarkony
,
J. Chem. Phys.
150
,
214101
(
2019
).
44.
X.
Zhu
and
D. R.
Yarkony
,
J. Chem. Phys.
137
,
22A511
(
2012
).
45.
Y.
Wang
,
C.
Xie
,
H.
Guo
, and
D. R.
Yarkony
,
J. Phys. Chem. A
123
,
5231
(
2019
).
46.
L. M.
Ibele
and
B. F. E.
Curchod
,
Phys. Chem. Chem. Phys.
22
,
15183
(
2020
).
47.
J. R.
Duke
and
N.
Ananth
,
J. Phys. Chem. Lett.
6
,
4219
(
2015
).
48.
F. A.
Shakib
and
P.
Huo
,
J. Phys. Chem. Lett.
8
,
3073
(
2017
).
49.
X.
Gao
,
M. A. C.
Saller
,
Y.
Liu
,
A.
Kelly
,
J. O.
Richardson
, and
E.
Geva
,
J. Chem. Theory Comput.
16
,
2883
(
2020
).
50.
S. J.
Cotton
,
K.
Igumenshchev
, and
W. H.
Miller
,
J. Chem. Phys.
141
,
084104
(
2014
).
51.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Theory Comput.
12
,
983
(
2016
).
52.
W. H.
Miller
and
S. J.
Cotton
,
Faraday Discuss.
195
,
9
(
2016
).
53.
H.
Meyera
and
W. H.
Miller
,
J. Chem. Phys.
70
,
3214
(
1979
).
54.
G.
Stock
and
M.
Thoss
,
Phys. Rev. Lett.
78
,
578
(
1997
).
55.
M.
Thoss
and
G.
Stock
,
Phys. Rev. A
59
,
64
(
1999
).
56.
H. D.
Meyer
and
W. H.
Miller
,
J. Chem. Phys.
71
,
2156
2169
(
1979
).
57.
U.
Müller
and
G.
Stock
,
J. Chem. Phys.
111
,
77
(
1999
).
58.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
145
,
144108
(
2016
).
59.
W. H.
Miller
and
S. J.
Cotton
,
J. Chem. Phys.
145
,
081102
(
2016
).
60.
N.
Bellonzi
,
A.
Jain
, and
J. E.
Subotnik
,
J. Chem. Phys.
144
,
154110
(
2016
).
61.
W. H.
Miller
and
S. J.
Cotton
,
J. Chem. Phys.
142
,
131103
(
2015
).
62.
G.
Tao
,
J. Phys. Chem. C
118
,
17299
(
2014
).
63.
S. J.
Cotton
,
R.
Liang
, and
W. H.
Miller
,
J. Chem. Phys.
147
,
064112
(
2017
).
64.
D.
Hu
,
Y.
Xie
,
J.
Peng
, and
Z.
Lan
,
J. Chem. Theory Comput.
17
,
3267
(
2021
).
65.
T.
Van Voorhis
,
T.
Kowalczyk
,
B.
Kaduk
,
L.-P.
Wang
,
C.-L.
Cheng
, and
Q.
Wu
,
Annu. Rev. Phys. Chem.
61
,
149
(
2010
).
66.
J. E.
Subotnik
,
E. C.
Alguire
,
Q.
Ou
,
B. R.
Landry
, and
S.
Fatehi
,
Acc. Chem. Res.
48
,
1340
(
2015
).
67.
A.
Kubas
,
F.
Hoffmann
,
A.
Heck
,
H.
Oberhofer
,
M.
Elstner
, and
J.
Blumberger
,
J. Chem. Phys.
140
,
104105
(
2014
).
68.
X.
Zeng
,
X.
Hu
, and
W.
Yang
,
J. Chem. Theory Comput.
8
,
4960
(
2012
).
69.
A.
Sirjoosingh
and
S.
Hammes-Schiffer
,
J. Chem. Theory Comput.
7
,
2831
(
2011
).
70.
F.
Webster
,
P. J.
Rossky
, and
R. A.
Friesner
,
Comput. Phys. Commun.
63
,
494
(
1991
).
71.
F.
Plasser
,
M.
Ruckenbauer
,
S.
Mai
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
,
J. Chem. Theory Comput.
12
,
1207
(
2016
).
72.
G.
Granucci
,
M.
Persico
, and
A.
Toniolo
,
J. Chem. Phys.
114
,
10608
(
2001
).
73.
P. O.
Löwdin
,
J. Chem. Phys.
18
,
365
(
1950
).
74.
L.
Joubert-Doriol
and
A. F.
Izmaylov
,
J. Chem. Phys.
148
,
114102
(
2018
).
75.
D. A.
Fedorov
and
B. G.
Levine
,
J. Chem. Phys.
150
,
054102
(
2019
).
76.
S.
Hammes-Schiffer
and
J. C.
Tully
,
J. Chem. Phys.
101
,
4657
(
1994
).
77.
N.
Yu
,
C. J.
Margulis
, and
D. F.
Coker
,
J. Phys. Chem. B
105
,
6728
(
2001
).
78.
L. G. C.
Rego
and
V. S.
Batista
,
J. Am. Chem. Soc.
125
,
7989
(
2003
).
79.
G. A.
Meek
and
B. G.
Levine
,
J. Chem. Phys.
145
,
184103
(
2016
).
80.
G. A.
Meek
and
B. G.
Levine
,
J. Chem. Phys.
144
,
184109
(
2016
).
81.
S.
Fernandez-Alberti
,
D. V.
Makhov
,
S.
Tretiak
, and
D. V.
Shalashilin
,
Phys. Chem. Chem. Phys.
18
,
10028
(
2016
).
82.
D. V.
Makhov
,
C.
Symonds
,
S.
Fernandez-Alberti
, and
D. V.
Shalashilin
,
Chem. Phys.
493
,
200
(
2017
).
83.
S.
Mai
,
M.
Richter
,
M.
Heindl
,
M. F. S. J.
Menger
,
A.
Atkins
,
M.
Ruckenbauer
,
F.
Plasser
,
L. M.
Ibele
,
S.
Kropf
,
M.
Oppel
,
P.
Marquetand
, and
L.
González
, SHARC2.1: Surface Hopping Including Arbitrary Couplings—Program Package for Non-Adiabatic Dynamics,
2019
, sharc-md.org.
84.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
, and
M.
Schütz
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
242
(
2012
).
85.
J. P.
Dahl
and
M.
Springborg
,
J. Chem. Phys.
88
,
4535
(
1988
).
86.
R.
Schinke
,
Photodissociation Dynamics: Spectroscopy and Fragmentation of Small Polyatomic Molecules
(
Cambridge University Press
,
1995
).
87.
D.
Tannor
,
Introduction to Quantum Mechanics: A Time-Dependent Perspective
(
University Science Books
,
Mill Valley
,
2007
).
88.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1370
(
2018
).
89.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Int. J. Quantum Chem.
115
,
1215
(
2015
).
90.
S. J.
Cotton
and
W. H.
Miller
,
J. Chem. Phys.
150
,
104101
(
2019
).
91.
P.
Vindel-Zandbergen
,
L. M.
Ibele
,
J.-K.
Ha
,
S. K.
Min
,
B. F. E.
Curchod
, and
N. T.
Maitra
,
J. Chem. Theory Comput.
17
,
3852
(
2021
).

Supplementary Material