Quantum dynamics of excited states involving processes such as charge or exciton transfer, [1][2][3] nonradiative recombination, 4,5 charge carrier transport [6][7][8] and trapping 9 as well as photoinduced chemical reactions is ubiquitous in nature and materials research. Such processes play a determining role in realizing molecular grounds of vision, 10 DNA photodamage protection, 11,12 photosynthesis 13 and artificial photocatalysis, 14,15 photovoltaics 16,17 and so forth. The nonadiabatic molecular dynamics (NA-MD) method [18][19][20][21] has been instrumental for simulating such quantum dynamical processes. While fully quantum NA-MD calculations are out of our reach due to their exponential complexity, except for lowdimensional problems or a select set of special problems, the approximate NA-MD schemes, particularly the family of the trajectory surface hopping (TSH) methods [22][23][24][25][26] have become a popular choice for modeling complex abstract and atomistic systems. Out of the family of the TSH schemes, the seminal Tully's fewest-switches surface hopping method (FSSH), 27 has become one of the most recognized and widely adopted approach due to the simplicity of its implementation as well as due to its improved ability to describe branching events and quantum-classical thermal equilibrium compared to traditional Ehrenfest dynamics. 28 Despite the general success of the traditional FSSH and alike TSH schemes, they still contain ad hoc features, defining some of their limitations, 23,29 such as the momentum jump approximation and subsequent momentum rescaling in the TSH. While the momentum rescaling in the direction of the nonadiabatic coupling vectors (NACVs) originating from the interpretation of Pechukas force 30 can be justified, [31][32][33] practical difficulties of obtaining NACVs in many situations stimulated the development of approximate, but truly ad hoc rescaling algorithms along the directions of nuclear momenta or state gradient differences. 34,35 Moreover, the variety of possible recipes for handling nuclear momenta in NA-MD simulations becomes even greater when frustrated hops are factored in. Frustrated surface hops arise when nuclear kinetic energy is insufficient to satisfy the total energy conservation when potential energy change upon the proposed state transition is too large. In this situation, nuclear momenta can be reversed, 36 kept unchanged, 37 or conditionally reversed. 38,39 Furthermore, although the TSH schemes are stochastic by definition, the discontinuity of the "classical" nuclear variables due to momentum rescaling appears an awkward concept, given that the starting equations are continuous. This philosophical difficulty stimulated a number of approaches where the need for discontinuous renormalization of dynamical variables is minimized to certain extent. [40][41][42][43][44][45] Another conceptual limitation of the FSSH-like schemes is the expectation that total energy conservation is satisfied at the level of individual trajectories. Considering that individual trajectories in quantum-classical methods may not always have physical meaning but should rather be treated as auxiliary objects, imposing energy conservation at the individual trajectory is an overly strict and not justified condition: what needs to be conserved is the ensemble average of the total energies among trajectories, not a single-trajectory energy. In this regard, an analogy to density functional theory (DFT) may be drawn. While the Kohn-Sham orbitals are often used in such calculations, they are not considered 1-electron states, but rather as the auxiliary objects needed to represent the overall charge density. Breaking the energy conservation at the individual trajectory level and allowing energy exchange between quantum-classical trajectories becomes especially fruitful for modeling scattering in superexchange process, where the nonadiabatic transition between two states A and B is mediated by another higher-energy state C. 46 If the states A and B are not directly coupled to each other but are coupled to the state C, the transition between them is prohibited when nuclear kinetic energy is low and all proposed hops are frustrated. By allowing energy exchange between trajectories, one can overcome the energy barrier and lead to successful transitions. To addressing this strict energy conservation problem, several NA-MD methods of utilizing multiple trajectories have been proposed, including second-quantized surface hopping (SQUASH) 47 and coupled trajectories mixed quantum-classical approach with energy-based decoherence (CTMQC-E) derived from the exact factorization (XF) framework. 48 To address the above conceptual flaws of the FSSH and FSSH-alike methods, quantum-trajectory surface hopping (QTSH) approach was developed by Martens and co-workers a while ago. 49,50 QTSH is derived from the QCLE 51,52 by applying the Wigner-Moyal transformation that maps operators to the corresponding phase-space functions leading to a systematic set of quantum-classical equations of motion. Unlike the TSH schemes, QTSH introduces an additional coherence energy beyond the diagonal energy from the active state, which naturally results in a quantum force that continuously adjusts the nuclear momenta before and after successful hops. Such forces act to conserve the total energy at the trajectory ensemble level yet allow the total energies of individual trajectories to vary. By construction, the momentum rescaling procedure is no longer needed. Despite the collective nature of the total energy conservation in QTSH, it is formulated as an independent-trajectory method, leading to its computational cost being comparable with that of the FSSH. It should be emphasized that while the QTSH solves the ad hoc momenta rescaling and trajectory-resolved total energy conservation problems, it still inherits the overcoherence problem known for FSSH. 33,53 This shortcoming was addressed in the early formulations of the QTSH of Martens by employing a simplified collapse-based decoherence correction.
To date, QTSH has been successfully applied to reproduce scattering probabilities in the Tully's single avoided crossing problem 27 as well as in the superexchange model. 46 Despite the conceptual attractiveness of the QTSH method, the prior works adhere to a two-state description. Furthermore, no practical implementation of the QTSH method in generalpurpose NA-MD software has been reported so far. Recently, Dupuy et al. presented an extension of the QTSH to multiple states and combining it with the XF surface hopping method (QTSH-XF). 54 However, the work utilized an ad hoc generalization of the QTSH theory, without a formal derivation of the corresponding equations or any generalpurpose software implementation.
To fill in the current gaps, we present a formal generalization of the QTSH method to multiple electronic states. Starting from the QCLE formalism, we derive the equations of motion for effective electronic and nuclear degrees of freedom (DOFs) which constitute the heart of the QTSH method. Our derivations lead to compact matrix-vector equations as well as the detailed analysis of several key terms in the equations (e.g., different contributions to forces), which enable assigning physical interpretation for such terms. In addition, our derivation suggests dissipative ("electronic damping") terms which have not been discussed before. We interpret them as the decoherence-causing terms and propose the combinations of the QTSH with decoherence correction algorithms based on simplified decay of mixing (SDM) 55 and XF, [56][57][58] denoted by QTSH-SDM and QTSH-XF. We implement the multistate QTSH approach in the open-source Libra package, [59][60][61] making it more accessible to the broader users' community. Employing the new implementation, we demonstrate that QTSH produces consistent results with those of FSSH, without relying on ad hoc momentum rescaling algorithms. We illustrate a linear correlation between energy conservation accuracy and internal consistency metric in QTSH. We find that applying decoherence corrections is not only needed to improve the internal consistency of simulations but is also essential for improving ensemble-level energy conservation within the QTSH algorithm.
In this section, we first review the core TSH algorithms and introduce the key notation. In most mixed quantum-classical approaches, nuclear density is represented by ensembles of classical trajectories. The corresponding nuclear degrees of freedom (DOFs) are characterized by nuclear coordinate q and momenta p which are propagated classically:
Here, the bold notation is for vector matrices containing the corresponding dynamical variables up to N dof , the number of DOFs:
T and so on. These vector-matrix notations are used throughout this paper. The dot over the symbols indicates the time derivatives.
Here, M is the diagonal matrix containing each mass of the DOF as its elements, i.e.,
), the index a appearing in eq 1b stands for the active state index of the trajectory considered. For brevity, we omit the trajectory index in this work with the understanding that all dynamical variables appearing in equations refer to a given trajectory, unless otherwise noted. The electrons are treated quantum mechanically-the electronic state is described by the electronic wave function Ψ(r, t; q), parametrically dependent on q, and evolving according to the time-dependent Schrodinger equation (TD-SE):
Here, ℏ is the reduced Planck constant and V ̂is the electronic Hamiltonian. In the practical implementation of the methods, we utilize the atomic unit system in which ℏ = 1. In most schemes, the time-evolving electronic wavefuction is expressed into a basis of adiabatic and stationary wave functions, {ψ i , i = 0,1, •••, N -1}, where N is the size of the adiabatic basis:
In the quantum-classical formulations, nuclei are treated classically, so the nuclear part of the electron-nuclear wave function is reduced to a time-dependent coefficients c i (t) for the corresponding electronic state ψ i . While the basis of stationary electronic states for TSH can be chosen rather flexibly, 62 the adiabatic representation has been argued as one of the most reliable. 63,64 In this work, we choose the adiabatic representation, {ψ i :V ̂ψi = V ii ψ i }, although the reformulation to the diabatic representation is straightforward.
Substituting eq 3 into eq 2, one obtains the corresponding equations driving evolution of electronic variables (amplitudes of the adiabatic basis functions in the time-dependent superposition):
where
is the time-derivative nonadiabatic coupling (tNAC) matrix element, and h ij = ⟨ψ i |∇ q ψ j ⟩ is the corresponding NAC vector (NACV). The electronic propagation, eq 4, can be written in terms of the density matrix, ρ ij = c i c j * or its real (α ij ) and imaginary (β ij ) components,
These equations describe the coherent evolution of the corresponding electronic subsystem. This evolution corresponds to a mean-field description of the process but can not account for quantum-mechanical branching effects, as well as it does not respect the detailed balance between electronic and nuclear subsystem, hence failing to correctly reproduce thermal equilibrium of quantum system.
2.2. Multistate QTSH Approach. While the typical TSH formulations start with the ad hoc separation of electronic and nuclear equations of motion, such equations can be derived more systematically starting from the quantum Liouville equation:
Here, the full system's Hamiltonian H ̂= T ̂+ V ̂acts on both electronic and nuclear coordinates, which are described by the density operator ρ. Here, T ̂is the nuclear kinetic energy operator, and V ̂is the electronic Hamiltonian operator also present in eq 2. We consider a complete basis orthonormal electronic states {|ψ i ⟩:
Substituting the resolution of identity condition to eq 6, one obtains a matrix form of this equation:
Here, H ̂ij = ⟨ψ i |H ̂|ψ j ⟩ and ρî j = ⟨ψ i |ρ|ψ j ⟩ are the matrix elements of H ̂and ρ̂in the electronic basis and are themselves operators belonging to the Hilber space of nuclear position states, |q⟩. In the Wigner-Moyal representation, 52,65 the nuclear operators becomes functions of the nuclear phase space coordinates, z = (q, p). Products of operators are represented by the star (or Moyal) product:
where q p p q = and the arrows indicate the direction in which the corresponding differential operators act. The quantum-classical Liouville equation (QCLE) is then obtained if one keeps only the lowest order nonclassical term in the expansion of the Moyal product:
where {A, B} = ∇ q A∇ p B-∇ p A∇ q B is the Poisson bracket of the phase-space functions A(z) and B(z). Applying the Wigner-Moyal transformation to eq 7, one obtains:
, ( , ) ( , ), ( ) ) ij k N ik kj ik kj k N ik kj ik kj 0 1 0 1 = + { } { } = = (10)
Here, all the quantities are now the phase-space functions:
(or as a superposition of auxiliary Gaussians), as shown in Sections S1 and S2 of Supporting Information, one arrives at the following equations for electronic and nuclear variables:
Here, F i = -∇ q V ii .
In the QTSH, the diagonal component of the force, eq 13 is reduced to an adiabatic (or diabatic) force of the active state
for a given trajectory, ∑ i F i ii →-∇ q V aa , assuming that the given trajectory k can at any time be only in one active state.
The derived equations contain numerically inconvenient gradients of the components of the NACV. The equations can be simplified by introducing the kinematic momentum P, proposed by Miller and co-workers, 50,66 which is defined as
Using eq 12, eq 14 becomes
Then the canonical and kinematic momenta have the following relationship.
, : , :
, :
Since
The time derivative of the kinematic momentum is extended further by inserting electronic equations of motion.
The equations of motion for the electronic variables are (see Section S4 of Supporting Information):
Here, ω ij = (V ii -V jj )/ℏ, g is an assumed spatial component of the density matrix function, ρ ij (z). Note that eqs 18a and 18b are nothing but eqs 5a-5c, just with added terms that can be regarded as "electronic damping" and cause decay of coherences. In other words, the extra terms on the RHS of eqs 18a and 18b can be interpreted as the terms causing electronic decoherence. The electronic damping term can be considered proportional to the rate of branching of nuclear basis functions (Section S9 of the Supporting Information). Assuming a Gaussian form of such functions, it may be related to the approach of Schwartz, Bitten, Prezhdo, and Rossky 67 which uses Gaussian overlap decay as a quantification of decoherence. In practice, we replace these terms with the ad hoc decoherence corrections of the two kinds, as discussed in Section 2.5.
Inserting eq 18b into eq 17, we have the force to evolve the kinematic momentum. Note that we have a new nonclassical force F (nc) as well as the usual active-state force F a in TSH.
, :
, : 0
Notably, the first-order force, F (1) , can be recognized as the off-diagonal components of the Ehrenfest force, and the form of the second-order force, F (2) , is akin to a Berry force, mediated by the density matrix. 68,69 The electronic friction 70,71 or decoherence force, F D , is led by the electronic damping in the electronic part.
One thing to keep in mind about eqs 19-23 is that the canonical momentum p is utilized in the equation relating the scalar NAC and NACV to each other in those derivations. However, replacing p with P in such derivations adds only a small contribution of the second order in ℏ, O(ℏ 2 ), and its effect is seen only modestly as discussed by Dupuy et al. 54 As further shown in Sections S5-S7 of Supporting Information, the nonclassical force can be written in terms of the vibronic Hamiltonian, so eqs 19-23 becomes
Here, we utilize the vector-matrix notation, where c = (c 0 , c 1 ,
In QTSH, the total energy of the electron-nuclear system, i.e., Tr(ρH) is represented through the phase-space densities:
Thus, total energy (eq 25a) contains the coherence energy (eq 25c) as well as the diagonal energy (eq 25b) from the stochastic hop. Consequently, the total energy can be simplified in terms of the kinematic momentum P k .
When the internal consistency condition, δ ̇i,ad k ≈ αi i,k , ∀ i, k, is satisfied, eq 26 yields E ̇tot = 0 if the evolution of electronic and nuclear DOFs is governed by the QTSH equations of motion discussed above (See Section S8 for the derivation of energy conservation in QTSH).
2.4. Hopping Algorithms. The coherent electronic evolution in TSH schemes, including QTSH, is replaced by the stochastic surface hops. The active state determined by such a procedure defines the nuclear force that drives the system. The hopping probabilities for FSSH and CSH (Consensus Surface Hopping) 51 are given:
Here, the σ function is defined as σ(x) = xH(x) = max (x,0) for the sake of conciseness. Here, H(x) is a Heaviside function. The CSH hopping probability scheme is similar to that of FSSH but introduces interaction of propagated trajectories. Here, the local densities ⟨ρ ii ⟩ n and coherence ⟨α ki ⟩ n are given as
As a result, the hopping probability for each trajectory depends on the state of all other trajectories, via the "consensus" among the entire ensemble. In QTSH, the ensemble average elements in eq 27b are replaced with "onsite" elements on each trajectory as QTSH is an independenttrajectory variation of CSH. This makes the hopping scheme in QTSH coincide with that from FSSH, eq 27a. We note that our earlier implementation of QTSH employed canonical momentum in the hopping probability in the adiabatic representation, rather than the kinematic momentum adopted here. This modification leads to more accurate agreement with exact quantum results. 50 One important point to note is that the capability of QTSH to capture detailed balance is yet to be fully investigated. While the FSSH hopping probability is known to achieve the detailed balance quite well, 72,73 its use within the QTSH framework may not guarantee the detailed balance. Additional studies of this matter would be desirable but go outside the scope of the current work.
In this work, we combine QTSH with the decoherence correction from the simplified decay of mixing (SDM) 55 and exact factorization (XF) [56][57][58] approaches, leading to QTSH-SDM and QTSH-XF. The direct use of the non-Hermitian electronic damping in eqs 18a and 18b, is out of scope of the current work, which could cause a numerical stability problem. Thus, we instead consider the decoherence effect with the pre-existing decoherence algorithms. The decoherence correction from SDM is characterized by the energy-based decoherence time τ ij . 74 i k j j j j j y
Here, |E i -E j | is the adiabatic energy gap for a given pair i and j, and E kin is the nuclear kinetic energy. The C parameter is an empirical parameter, and its default is 0.1 Ha. 75 According to the decoherence time, eq 29, coefficients for the nonactive state i ≠ a are damped by an exponential factor, exp(-Δt/τ ia ). Afterward, the coefficient for the active state is renormalized so that the norm is conserved. The decoherence correction from XF replaces the electronic damping term in the electronic propagation in eqs 18a and 18b. Then the resulting electronic equations become:
The decoherence terms are governed by the following XF Hamiltonian matrix:
Here, the quantum momentum and the phase gradient matrix ϕ ν are determined by the auxiliary trajectory on each adiabatic state i, (q i , p i ) as follows.
Here, σ ν in eq 32a is the width parameters for auxiliary wave packets and t i is the time when the auxiliary trajectory on state i is spawned. To determine the value of the σ ν parameter, one can adopt a characteristic fixed width or time-dependent width such as Schwartz 76,77 or Subotnik width. 78,79 For the details of the auxiliary trajectory propagation and the time-dependent width approximation, refer to ref 57. Notably, the XF decoherence in eq 30 has a similar form of the electronic damping, i t ( )
in eq 11 if g is assumed to be the phase-space density in the form of a Gaussian function, except that the electronic damping yields an additional momentum counterpart.
When SDM and XF decoherence corrections are applied to QTSH, only electronic propagation is modified in both cases. In QTSH-SDM, the damping and renormalization of coefficients are conducted after the electronic propagation. In QTSH-XF, the half-time evolutions for the XF Hamiltonian matrix, eq 31, are applied to electronic propagation, leading to symmetric trotterization of the whole electronic Hamiltonian, V + H XF . In this work, we neglect the decoherence force. Improving internal consistency is a major factor in energy conservation. Also, in QTSH-XF, the XF quantum force, characterized by the inverse mass factor, is known to have a minimal impact on the dynamics. 54 2.6. QTSH Implementation in the Libra Package. The multistate QTSH formalism discussed above is implemented in the Libra package. 59,60 While the key equations have been discussed already, we put several pragmatic comments here. First of all, the state tracking and phase consistency corrections of adiabatic states and the dependent properties are of critical importance for the formally well-defined equations to work as expected. In practice, trivial crossings are not uncommon, especially in large systems with many states. Here, the state crossing may result in an abrupt sign change of NACs and NACVs during a given integration time-interval. This problem can be solved by either utilizing the min-cost (aka Munkres-Kuhn or Hungarian) algorithm 80 in combination with the phase correction algorithm. 81 Alternatively, the local diabatization scheme 82,83 has shown great use in many previous applications. Both approaches find a basis reprojection matrix,
T, which takes care of both state permutations (trivial crossings) and sign/phase changes:
Here, tilde notation refers to the reordered and phasecorrected basis functions, while the regular notation indicates the "raw" (as obtained from the eigensolver procedure) adiabatic states. Also, the notation |ψ⟩ = (|ψ 0 ⟩, |ψ 1 ⟩,
•••, |ψ N-1 ⟩) is used. The computed reprojection matrix is then used to determine the active state permutation (to reassign the active state correctly), to transform the state coefficients (to preserve the continuity of the overall time-dependent wave function), c(t + Δt) = T -1 c(t + Δt). It should also be used to transform the electronic structure properties that are determined by the basis. For instance, if the following "raw" value of the matrix representation, A(t + Δt) = ⟨ψ(t + Δt)|A ̂| ψ(t + Δt)⟩ of an operator A ̂is known, it should be transformed to the reprojected basis |ψ̃(t + Δt)⟩ as Ù Ù Ù A A A t t t t t t A t t T t t A t t T T t t T ( )
Examples of such operators could be the derivative coupling tensor and Hamiltonian matrix. Considering this precaution, QTSH is implemented based on already existing general TSH workflow, but with the following changes:
• Electronic damping in eqs 18a and 18b is either neglected (in regular QTSH) or its effect is introduced in an ad hoc way by augmenting the QTSH prescription with the effective decoherence correction. As options, we consider the SDM and XF decoherence corrections, leading to QTSH-SDM and QTSH-XF methods, respectively. Adding the dissipation (electronic "damping") term in the RHS of eqs 18a and 18b directly will be addressed in future works; • All proposed hops are now accepted. Hence, there are no frustrated hops, and there is no ambiguity regarding what to do if such hops occurred (as in the original FSSH approach). • No momentum rescaling is utilized upon successful hops. Energy conservation is now expected to be satisfied only at the trajectory ensemble level, not for individual trajectories and is maintained due to the presence of nonclassical forces. Since our current implementation replaces the last terms on the RHS or eq 11, or, equivalently, eqs 18a and 18b, with the ad hoc decoherence corrections, the total energy conservation is satisfied only approximately, even in the ensemble sense. • Nuclear dynamics is now affected by the nonclassical forces except for the decoherence force, eqs 19-23; • The nuclear dynamics is conducted in terms of kinematic momenta-the canonical ones are computed from the kinematic, if necessary. • The hopping probability in the current QTSH implementation is identical to the hopping probability of FSSH, eq 27a.
Hamiltonians and Computational Details. Holstein, superexchange, and phenol models (Figure 1) are chosen to evaluate the performance of the family of QTSH methods, i.e., QTSH, QTSH-SDM and QTSH-XF compared to that of FSSH and to validate our implementation of these methods in the Libra package. 59,60 The 1D Holstein model 84 is a bound potential with two diabatic surface crossings. It highlights the dynamics with multiple crossings, which can lead to nontrivial branching, decoherence, and interference phenomena. The diabatic Hamiltonian for 1D Holstein model is defined as
with E 0 = 0, E 1 = -0.01, R 0 = 0, R 1 = 0.5, k 0 = 0.002, k 1 = 0.008, and V = 0.001. All parameters are given in the atomic unit. The superexchange model 46 is chosen as a wellestablished 1D multistate example with several references available. 46,49,85 We also utilize the phenol model Hamiltonian for the photoinduced hydrogen elimination reaction 86 for exploring a 2D potential landscape, described by the OH length r and the CCOH dihedral angle θ. The parametrization for the superexchange and phenol models can be found in the corresponding references.
For benchmarking. the quantum dynamics is performed with the DVR method. 87 The initial wavepacket χ(q,0) is set to the following Gaussian wavepacket. i k j j j j j j j y { z z z z z z z q s q q s ip q q ( , 0)
)
Here, q 0,f and p 0,f = ℏk 0, f are the average position and momentum for the
), with the corresponding Gaussian width, s q,f . The standard deviations of q and p, corresponding to the probability density |χ(q,0)| 2 are given as s / 2 q f q f , , = and σ p,f = ℏ/2σ q,f . For the 1D Holstein model, the mass is 2000 au. The initial wavepacket is on the ground state with q 0 = -4.0 Bohr at rest (p 0 = 0 au) and σ q = 1.0 Bohr. The DVR simulation is conducted for 16,000 au with the time step of Δt = 1.0 au and a uniform grid width of 0.025 Bohr in range [-25, 26.175] Bohr (2048 grid points). For the 1D superexchange model, the mass is 2000 au The initial wavepacket is on the ground state with q 0 = -10.0 Bohr and σ q = 1.0 Bohr. The initial momentum ℏk 0 varies in the interval [3.0,20.0] au for the scattering calculations. The DVR simulation is conducted for 40,000 au with the time step of Δt = 1.0 au and a uniform grid with the coordinate grid point spacing of 0.025 Bohr and the grid range to be [-150, 259.575] Bohr (total of 16,384 grid points). For the 2D phenol model, the reduced mass for r is 1728.5 au, and moment of inertia for θ is 5132.0 au The initial wavepacket is on the first excited state from the vertical transition at the equilibrium position on the ground state potential energy surface (PES), defined with q 0 = (r 0 , θ 0 ) = (0.96944 A°,0 rad), p r,0 = 15 au, and ( , ) (0.092/ 2 Å, 0.55/ 2 rad)
, following the setting of the work of Pollien et al. 86 The DVR calculation is conducted for 4,000 au with Δt = 10.0 au and the grid widths of Δr = 0.02 Bohr and Δθ = 0.03 rad in range [0, 81.9] Bohr and [-60, 62.85] rad, respectively (4096 grid points for each). The grid for θ is deliberately extended instead of using the periodic boundary conditions. In all DVR calculations, the initial state is defined in the adiabatic representation, consistently with the initialization of the TSH calculations. Since the DVR calculations are conducted in the diabatic representation, the diabatic-to-adiabatic transformation is used to project the initial adiabatic wavepacket on the corresponding wavepackets (potentially on multiple states) in the diabatic representation. The split operator Fourier transform integration is conducted in the diabatic representation, and the corresponding adiabatic properties are computed according to the underlying transformation.
In all TSH calculations, the initial coordinates and momenta are sampled from the Gaussian distribution based on the wavepacket probability density for the DVR dynamics: i k j j j j j j j y { z z z z z z z i k j j j j j j j y { z z z z z z z q p q p P q q p p ( , ; , ,
( )
For FSSH, the Jasper-Truhlar criterion 38 is applied when the dynamics encounter frustrated hops. 2000 trajectories are used to obtain statistical results.
The comparison between the TSH and DVR dynamics is made in terms of population, coherence indicator and the branching ratio. The total wave function Ψ(r, q, t) is expressed in the electronic basis {Φ i (r; q)} as r q q r q t t ( , , ) ( , ) ( ; )
where the total nuclear density |χ| 2 is given as the sum of densities of each adiabatic wave packet,
Then the population and coherence indicator in DVR are computed by integrations over the grid.
Interpreting the coefficient c i as c i = χ i /χ and approximating the total nuclear density to the summation of delta function centered at each trajectory, 88,89 the relevant population and coherence indicator expressions in the TSH dynamics become
In eq 40b, we used the property
The population from eq 40a is the so-called SE population, reflecting on the fact that the population is based on the electronic coefficients propagated by the Schrodinger equation, while the population from the active state is called the SH population.
P N
We compute both population expressions to check the internal consistency between them.
For the comparison of the nuclear distribution between the TSH and DVR calculations, the transmission and reflection for [3.0, 20.0] au are computed for the superexchange model. To compute the reflection and transmission on each state with respect to the initial momenta, active states are counted in the TSH dynamics by dividing the grid space into (-∞, 0) and (0, ∞) at the end of simulation. In the DVR calculation, the branching ratio calculation is calculated by integrating each adiabatic wavepacket over space using the same divisions. For the phenol model, the dissociation probability based on the cutoff OH length r cutoff = 2.6 Å is computed by counting the trajectories where r > r cutoff from the TSH calculations and integrating the nuclear density with the limits of (r cutoff , ∞) and (-∞, ∞) from the DVR calculation.
To assess the methodologies in the case of multiple state crossings, we compute the dynamics in the Holstein model that can support multiple branching at double diabatic crossing regions (Figure 2). We observe that regular QTSH scheme causes systematic drift of the total energy due to accumulation of errors at every state crossing region (large nonadiabatic coupling regions), as shown in Figure 2a. The decoherencecorrected QTSH-SDM and QTSH-XF schemes notably reduce the total energy drift (Figure 2a) and lead to better energy conservation quality. The decoherence schemes also improve the internal consistency, making the SE and SH population agree better with each other (Figure 2b). In this regard, the QTSH-XF scheme yields the population dynamics in remarkable agreement with that of the DVR reference (Figure 2b), while the QTSH-SDM, QTSH, and FSSH methods all yield more distinct dynamics but comparable to each other.
The system repetitively passes the strong coupling regions, leading to large Rabi-like oscillations of the state populations, whose amplitude decays due to quantum mechanical branching and equilibration. As far as coherences are concerned, both QTSH-SDM and QTSH-XF schemes show the evolution of the coherence indicators in qualitative (and pretty close quantitative) agreement with that of the DVR calculations (Figure 2c). The bare QTSH scheme yields coherence indicators comparable to those of the overcoherent FSSH method (Figure 2c). This is expected since QTSH is not meant to introduce any decoherence effects on its own. The overcoherence of the bare QTSH scheme causes imbalance of the diagonal, eq 25b, and coherence, eq 25c, energies, resulting in the poor total energy conservation shown in Figure 2a. Incorporating effective decoherence corrections is crucial to achieve fair energy conservation. The correlation between energy conservation and internal consistency becomes more evident in the following section.
Model. First, we conduct the scattering calculation for the 1D superexchange model for a range of average initial nuclear momenta using the FSSH, QTSH, and decoherence-corrected QTSH schemes and compare them to the numerically exact results. For QTSH-SDM, the energy decoherence parameter is set to its default value of C = 0.1 Ha, while for QTSH-XF, the width parameter is defined as σ = 0.1 Bohr, corresponding to one-tenth of the width of the initial Gaussian distribution, consistent with the previous work from Dupuy et al. 54 The initial momentum range, spanning from 5.0 to 20.0 au, has been explored in the original QTSH work of Martens. 49 The low limit of this regime, ℏk 0 = 5.0 au, particularly important as it represents a above which no substantial reflection occurs during the scattering process (Figure 3). This is attributed to the PES of the superexchange model (Figure 1b), where the energetic barrier between the ground and the first excited state, 0.005 Ha, is lower than the kinetic energy at this threshold, (ℏk 0 ) 2 /2m = 0.00625 Ha. Consequently, classically forbidden hops between the ground and the first excited states are absent in this regime.
The transmission probabilities on states 0 and 1 for the high k regime (Figure 3a) computed using either the QTSH or its decoherence-corrected versions (QTSH-SDM or QTSH-XF) agree with those obtained from the FSSH and numerically exact calculations. This is a remarkable result considering that the velocity rescaling procedure is eliminated in QTSH by construction. This agreement serves as one of the validations of correctness of the method's implementation. In this domain, where the initial nuclear momentum exceeds 5.0 au, the majority of wavepacket prepared initially on the ground state traverses the coupling region only once (see the Supplementary movie 1). As a result, the overcoherent feature of electronic propagation from FSSH and QTSH does not affect the prediction of adiabatic populations significantly and therefore agrees well with the fully quantum calculations. Obviously, the regime of greater interest lies at lower momenta, where reflection contributes to the dynamics to a certain extent (Figure 3b,c), and where forbidden hops would be occurred in the original FSSH scheme, while the effect of natural wavepacket broadening 41,90 is minimal. This regime is explored in the following sections.
3.2.2. Energy Conservation in the Low k Regime in QTSH. We now focus on the QTSH dynamics in the low k regime, ranging from 3.0 to 5.0 au still using the superexchange model. Under such conditions, reflection becomes apparent, and the nuclear wavepacket gains an additional channel to propagate (Figure 3c). The wavepacket, initially on the ground state, branches into the ground and the first excited states at the coupling region. The portion of the wavepacket remaining on the ground state continues to traverse the ground-state PES. Meanwhile, the wavepacket on the first excited state is temporarily trapped and undergoes multiple nonadiabatic transitions into the ground state (See the Supplementary movie 2).
When both reflection and transmission occur during the dynamics such as in the case of ℏk 0 = 4.0 au (Figure 4), the internal consistency of the TSH method breaks down, meaning that SE-and SH-based state populations become distinct (e.g., Figure 4c). This is in striking contrast to the case of ℏk 0 = 5.0 au (Figure 4d), where only transmissions occurred, and no internal consistency was broken. Expectedly, for both initial momenta, ℏk 0 = 4.0 and ℏk 0 = 5.0, the overcoherence of the QTSH is manifested via nondecaying values of the coherence indicators, in contrast to exact quantum simulations (Figure 4e,f). Notably, with the lower initial momentum, overcoherence becomes more drastic because the trajectories on the first excited state encounter the coupling regions multiple times, developing stronger coherences between the ground and first excited states. On the other hand, the SH populations, which depend on the discrete active state, are more reliable in such scenario and yield a closer agreement with the DVR reference (Figure 4c).
Together with the worsening of the internal consistency for lower initial momenta, we observe more notable violation of the energy conservation (Figure 4a). Compared to the case of ℏk 0 = 5.0 au (Figure 4b), where the balance between the diagonal energy and coherence energy is well maintained, for the ℏk 0 = 4.0 au it is broken. This imbalance appears when the SE and SH populations start to deviate, leading to a drop in the total energy. Eventually, the total energy is stabilized after all trajectories move away from the coupling region, since the coherence energy is governed by NAC, eq 25c, and hence vanishes in the asymptotic regions. This observation demonstrates that the energy conservation at the trajectory ensemble level is affected by the quality of the internal consistency of the SE and SH population as discussed in Section 2.3.
To investigate the relationship between the quality of the internal consistency of SE and SH populations and the total energy conservation, we calculate root-mean-square errors (RMSEs) of the SE and SH population difference (eqs 40a and 41) and total energy fluctuations in the lower-k regime (Figure 5). A clear linear correlation is observed between such quantities. The total energy fluctuation increases as the initial nuclear momenta decrease, aggravating the overcoherence due to smaller tNAC. Thus, addressing the internal consistency problem is crucial for obtaining reliable results using QTSH.
Our rationalization of the observed correlation is the following. The rigorous derivation of the QTSH equations of motion yields the "dissipation/electronic damping" term where Σ is the diagonal matrix containing the Gaussian width for each phase-space DOF, and Δz k is the deviation of phase space variables, z k -z ̅ (see Sections S9 of Supporting Information). Thus, the "electronic damping" term is closely related to branching of the nuclear dynamics. In the high k regime, the transmission is predominant and there is no significant difference between individual trajectory's z k values and average z ̅ values, leading to small Δz k and hence to small dissipation terms. Thus, the effect of "electronic damping" in this domain is negligible. On the other hand, in the low k regime, reflection also contributes to the dynamics, and the phase-space deviation Δz k increases along with the branching. Thus, in the low k regime, correct treatment of the decoherence effect becomes a more crucial factor. 3.2.3. Combination of QTSH and Decoherence Corrections. As discussed in the previous section, one needs to introduce a decoherence correction into QTSH to improve the internal consistency and hence improve total energy conservation. In this work, we consider a set of widely used decoherence corrections instead of using the "dissipation/ electronic damping" terms i t ( )
in, eq 11, directly.
The consideration of such term in equations of motion leads to nonunitary dynamics, which may be more difficult to analyze and goes outside the scope of this work. In addition, since the nuclear envelope, g, is arbitrary, a nonunique choice of such terms may be present. By construction, the dissipation term can dampen coherences and hence shall be responsible for capturing decoherence effects. Thus, instead of using the term directly in eq 11, we combine the QTSH formulation with this term in eq 11 omitted with some of the well-known decoherence corrections, such as SDM and XF, denoted as QTSH-SDM and QTSH-XF, respectively. It turns out that in the lower-k regime, a careful choice of decoherence parameters is necessary to improve energy conservation as well as the internal consistency for each value of ℏk 0 , rather than using a uniformly defined parameters as in cases with larger initial momenta (Section 3.1).
We observe that the QTSH-SDM method improves the quality of the total energy conservation in the superexchange model dynamics at ℏk 0 = 4.0 au, when the energy decoherence parameter is set to C = 0.01 Ha (Figure 6a). Notably, the internal consistency of SE and SH populations is a necessary but not sufficient condition for ensuring total energy conservation, as shown by QTSH-SDM calculations with C = 0.01 Ha and C = 0.1 Ha (Figure 6a,b, respectively). As the C parameter is increased, the decoherence time also increases, and the dynamics becomes more coherent, leading to poorer energy conservation. However, the decoherence corrections such as the instantaneous decoherence approximation (IDA), 91 where the wave function is collapsed to a pure state at every attempted state hop, are not necessarily superior. When IDA approach is combined with QTSH, it eliminates the coherence energy budget (eq 25c) exclusively for trajectories undergoing hops, while keeping the overcoherence for trajectories that have not attempted surface hops. This imbalance eventually disrupts the whole energy compensation in QTSH (Figure S2). Unlike the IDA, the removal of overcoherences by the SDM and XF schemes need not be informed by the occurrence of surface hopping events, and hence the latter treat all trajectories on the same footing.
Selecting a reasonable SDM C parameter is also important to capture the correct dynamics of coherence near NAC region while also maintaining the internal consistency. Comparing the 6a,c). However, precautions need to be taken, as the coherence computed by the TSH method could contain intrinsic inaccuracies due to the trajectory-based approximation. The DVR values serve only as a qualitative guide for choosing the parameter. The "optimal" decoherence parameter, which achieves both energy conservation and internal consistency, should be carefully examined along with trajectory convergence (Figure S1).
We also present the QTSH-XF results for various width parameters when ℏk 0 = 5.0 au (Figure 7). As the XF decoherence correction involves ballistic motion of the auxiliary trajectories spawned on all trajectories, 57 for values ℏk 0 < 5 this motion easily becomes unstable/unphysical. Hence, for QTSH-XF approach, we analyze the dynamics with the initial conditions different from those considered for QTSH-SDM. In this case, energy conservation is less sensitive to the width of the auxiliary wavepackets σ, which controls the decoherence rates (Figures 7a-c and S3). As in the previous QTSH calculations (Figure 4b,d,f), when ℏk 0 ≥ 5.0 au, the trajectories pass through the coupling region only once, characterized by the flattened SH populations after a significant population exchange at t = ∼100 fs (Figure 7d-f). This indicates that subsequent inaccuracies of internal consistency do not impact the total energy conservation, as the nuclei would be already outside the coupling region, and the coherence energy, eq 25c, would become zero. Only early coherence around t = ∼100 fs has a finite contribution to the coherence energy.
3.3. 2D System: Phenol Hamiltonian. We now investigate the nonadiabatic dynamics of the phenol model for the photoinduced hydrogen elimination, using QTSH, QTSH-SDM, and QTSH-XF methods to evaluate their performance for a 2D problem. The nuclear wavepacket initially on the first excited state may transfer to the second excited state and remain transiently bound there or it can continue evolving along the first excited and ground state PES (Figure 1c) leading to molecular dissociation.
Populations from all TSH methods follow the overall DVR trend, despite the deviation from the reference for longer evolution times (Figure 8c,e,g). FSSH and QTSH exhibit comparable behavior, both failing to account for internal consistency. This result is similar to the one discussed for the superexchange model but exaggerated by multiple crossings of the strong nonadiabatic coupling regions. The SH populations produced by FSSH and QTSH methods are in agreement with each other. Similar to the 1D problem, the decoherence corrections improve the internal consistency of the methods, making SH and SE populations behave in a similar way (Figure 8c,e,g). The enhanced internal consistency from QTSH-SDM and QTSH-XF also leads to better energy conservation than for original QTSH, although the QTSH-XF method also leads to the QTSH-level total energy drift for longer simulation times (Figure 8a). Considering the similarity of the SH populations among the TSH methods, utilizing the SH populations is a preferred way to obtain more reliable descriptors that partially bypass the overcoherence problem. 92 In all TSH calculations, the population of the second excited state decays to a larger asymptotic value compared to the DVR reference, leading to the underestimation of the dissociation probability (Figure 8b).
Notably, the coherence indicators (Figure 8d,f,h) provide qualitative explanations for the SE populations. For all TSH methods, as coherence deviates from the DVR results at the early stage of the dynamics, overall populations start to deviate considerably from the DVR results at t = ∼15 fs. FSSH and QTSH show overcoherence, while the QTSH-SDM and QTSH-XF yields the "undercoherence". Notably, the coherences related to the ground state, ⟨|ρ 02 | 2 ⟩ and ⟨|ρ 01 | 2 ⟩, explains better alignment of the ground-state SE populations from FSSH and QTSH with the DVR reference compared to those from QTSH-SDM and QTSH-XF. The individual behavior of the coherences related to the ground state, that is, ⟨|ρ 02 | 2 ⟩ and ⟨|ρ 01 | 2 ⟩, deviates from the exact one. However, their summation is close to the DVR reference, leading to a cancellation of errors on average. For QTSH-SDM and QTSH-XF, coherences are underestimated starting from t = ∼20 fs (Figure 8g,h), and the degree of the ground-state relaxation starts to be overestimated compared to the DVR reference from that point.
In this work, we present the QTSH formalism generalized to an arbitrary number of states, including its detailed derivation, as well as its first implementation in the open-source Libra package. The implementation is validated through modeling nonadiabatic dynamics in 1D Holstein, superexchange, and 2D phenol models. As anticipated, the QTSH method produces results consistent with those of the FSSH scheme, although without the need for the ad hoc velocity rescaling procedure and hence without frustrated hops present in FSSH. We also combine the QTSH algorithm with the SDM and XF decoherence corrections, resulting in the QTSH-SDM and QTSH-XF, respectively. We show that using decoherence correction with the QTSH is important for ensuring its internal consistency and significantly improving the quality of the total energy conservation. Our calculations suggest a strong linear correlation between the degree of internal consistency (as measured by the SE and SH population difference) and the quality of the total energy conservation (as measured by the relative total energy fluctuation).
The QTSH scheme and especially its decoherence-corrected QTSH-SDM and QTSH-XF schemes show a reasonable agreement with the reference calculations for the current set of model Hamiltonians. The bare QTSH generally fails to conserve the total energy, even in the trajectory ensemble sense, when the branching effects are notable such as in the Holstein, phenol, and low-momentum superexchange models. This effect originates due to the current neglect of electronic "damping" terms and can be mitigated by including ad hoc decoherence corrections. Out of the two decoherence correction schemes combined with the QTSH, we find that QTSH-SDM is generally more robust as far as the total energy conservation is concerned, but QTSH-XF often yields a better agreement of population and coherence indicator dynamics with those from fully quantum calculations. These results highlight the need for ongoing investigation of the general problem of quantum coherence in mixed quantum-classical methods. The energy conservation and population/coherence predicted by the QTSH-SDM and QTSH-XF methods can be further controlled by corresponding hyperparameters of the methods.
The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.4c01751.
The PDF document containing derivations of the QTSH formalism, trajectory convergence tests, QTSH-IDA dynamics in the superexchange model, total energy fluctuation in the superexchange model with ℏk 0 = 5 au; DVR dynamics movies for the superexchange model with ℏk 0 = 5 and 4 au The Python scripts for all calculations and Jupyter notebook for data visualization are available via Zenodo (10.5281/zenodo.14538540). The QTSH algorithm in this work is available in the latest GitHub repository of Libra (https://github.com/ Quantum-Dynamics-Hub/libra-code/tree/devel commit 3adb48fb8732a6190656e9e118af7f9961931534) (PDF)
Wavepacket dynamics for the initial nuclear momentum of 5.0 au, with the majority in the ground state crossing the coupling region once (MP4) Wavepacket dynamics for the initial nuclear momentum of 4.0 au, with the wavepacket temporarily trapped on the first excited state and undergoing multiple nonadiabatic transitions to the ground state (MP4)
https://doi.org/10.1021/acs.jctc.4c01751 J. Chem. Theory Comput. 2025, 21, 2839-2853
Journal of Chemical Theory and Computation