Energy-conserving coupled trajectory mixed quantum–classical dynamics

The first law of thermodynamics states that the energy of a closed system must be conserved. For molecular systems, this means energy exchange must occur between the nuclei and electrons in such a way that the total energy of the molecule is invariant. To simulate the dynamics of complex molecules, approximations are inevitably required, and when energy conservation does not arise naturally in a given approximation, it is usually imposed as an extra condition. In some cases, the energy non-conservation is a consequence of simplifications made in the numerical implementation of the method, e.g., in the independent-trajectory version of the multiconfigurational Ehrenfest method in finite-Gaussian basis set implementations.^1,2 Justified by the large nuclear–electron mass ratio, mixed quantum–classical (MQC) schemes are the basis for many approximations for coupled electron–nuclear dynamics, where one propagates an ensemble of classical nuclear trajectories, each associated with a set of quantum electronic coefficients of, typically, a Born–Oppenheimer (BO) basis. Ehrenfest dynamics (EH) and surface-hopping (SH)^3,4 are the most widely used MQC methods, with SH generally preferred due to its ability to capture wave-packet splitting after regions of coupling between BO surfaces are encountered. Both conserve energy, albeit in different ways. Both involve a nuclear force that is the gradient of an electronic expectation value, with the electronic equations having a Hamiltonian evolution. However, while the consistency in the electronic and nuclear equations in Ehrenfest yields energy conservation on an individual trajectory level, the inconsistency in having coherent electronic evolution while each nuclear trajectory is collapsed on a given surface at any time means that energy conservation in SH is a more subtle issue.

SH imposes energy conservation on each individual trajectory by rescaling the velocities after a hop has occurred such that the gain or loss in potential energy is compensated by a loss or gain in kinetic energy. Not only is there no unique way to do this,^5–7 but the notion itself lacks proper physical justification since it is the energy of the ensemble of trajectories that mimics the underlying quantum wave packet that should be conserved, not the energy of an individual trajectory, and this too tight constraint leads to inaccuracies. One example is the so-called frustrated hop, where strict energy conservation precludes dynamical pathways that are accessible in the full quantum mechanical treatment. Furthermore, frustrated hops exacerbate the difference between electronic populations predicted from the electronic coefficient evolution and those predicted by the fraction of nuclear trajectories on the corresponding electronic state, which is known as the internal consistency problem of SH. More sophisticated SH schemes that do conserve energy over the ensemble rather than for individual trajectories^8,9 have not been applied beyond model systems so far. However, even aside from the incorrect imposition of individual trajectory energy conservation and its ad hoc manner of imposition, there is a third issue: if the electronic populations in the definition of the energy were obtained from the coherent electronic evolution instead of the fraction of trajectories, then SH methods would, in fact, not conserve total energy. Although SH schemes that include a decoherence correction^10–13 ameliorate the internal inconsistency (albeit in an ad hoc way), the other issues remain.

The exact factorization approach (XF)^14–23 paved the way for the development of new MQC methods that tackle the issues of standard MQC schemes. The coupled-trajectory mixed quantum–classical algorithm (CTMQC)^24–26 was derived by taking the classical limit of the XF equations and approximating some of the terms in ways justified by exact studies.²⁷ Electronic decoherence and nuclear wave-packet splitting emerge naturally, and CTMQC has successfully been applied to a number of interesting photo-induced dynamics.^26,28–32 The XF terms introduce new mechanisms of population transfer driven by nuclear quantum momentum, shown to be crucial in capturing dynamics of multi-state intersections²³ in an algorithm where the electronic equation was used in a SH scheme.^23,33–39

However, due to the approximations made to the XF equations, the ensemble of coupled trajectories in CTMQC is not guaranteed to conserve the total energy. In this work, we analyze why, and propose and demonstrate a modification, CTMQC-E, that restores energy conservation.

Our first example is the one-dimensional Tully’s extended coupling region (ECR) model,³ shown in the inset of Fig. 1. We study two situations where a Gaussian nuclear wave packet centered at −15 a.u. is incoming on the lower surface from the left with a high (k₀ = 32 a.u.) and a low (k₀ = 26 a.u.) initial momentum; the variance of the Gaussian wave packet at time zero is 20 times the inverse of the initial momentum. These two scenarios were studied with independent-trajectory XF-based methods in Ref. 41. For the higher initial momentum case after passing through the NAC region where some population is transferred to the upper surface, the wave-packet component on the lower surface moves faster separating in nuclear space from the part that has transferred to the upper state and the branches decohere. For the lower initial momentum, much of the wave packet transmitted to the upper surface reflect back and recross the NAC region, leading to a second splitting event. In our MQC simulations, 1000 Wigner-distributed trajectories with a fixed initial momentum k₀ were run starting on the lower surface, and the time-step used in the calculations was 0.1 a.u. Figures 1 and 2 show the excited state population $∑αρ22(α)$⁠, the fraction of trajectories on the active state for SH N₂(t)/N_tr, the coherence $∑α|ρ12(α)|2$⁠, and the energy from Eq. (10) for k = 32 a.u. and k = 26 a.u., respectively.

Figure 1 shows that for k₀ = 32 a.u. momentum, all MQC methods closely reproduce the exact populations. CTMQC and CTMQC-E give identical coherences that decay faster than the exact, whereas SH and EH do not decohere. The top panel shows that CTMQC-E conserves the total classical energy, curing the increase in CTMQC between 500 and 1000 a.u., where the quantum momentum is active. SH and EH conserve energy so are not shown in the top panel. After 1000 a.u., the trajectories fully decohere and energy, thereafter, remains constant in both methods. These results could be compared with the independent trajectory XF-based schemes shown in Fig. 5 of Ref. 41. There, the populations and coherences are also closely reproduced, but there are some variations, depending on the way the width parameter of the auxiliary trajectory is computed, and the equations are integrated.

Figure 2 shows the k₀ = 26 a.u. case. Consider first the standard MQC schemes. EH does not capture any reflection or decoherence after the first passage through the NAC region. On the other hand, the fraction of trajectories population measure of SH predicts an almost complete reflection of all trajectories on the upper surface distinct from the behavior of its electronic populations, displaying a poor internal consistency. As discussed earlier, if we were to take the electronic populations as the population observable, SH would not conserve energy. Now, let us turn to the XF-based MQC methods. CTMQC-E gives an improvement over CTMQC in the populations, coherence, and total energy. We note that to avoid numerical instabilities occurring when the velocity becomes too small in the denominator of Eq. (14) (an issue for reflecting trajectories), we apply a velocity threshold such that the new definition of accumulated force from Eq. (14) is used only above a small threshold (10⁻⁵), and otherwise, the original definition Eq. (7) is used. This procedure can lead to small discontinuities in the energy, as seen in the red curve in Fig. 2 at around 2000 a.u. when the energy jumps up a little above the initial value. Further work will involve stabilizing the algorithm, but we see already the improvement in energy conservation in CTMQC-E over CTMQC. Furthermore, the populations are slightly improved and also the coherence. These results may again be compared with the independent trajectory XF-based schemes in Fig. 4 of Ref. 41, where the results were quite sensitive to details in the algorithm, and a time-dependent auxiliary-trajectory width was needed to capture the exact dynamics in the energy-conserving scheme.

These three modes, which drive the ultrafast cis–trans isomerization of rPSB11,^47–50 also drive the isomerization of cis-PSB3 after photo-excitation from the ground (S₀) to the first singlet excited state (S₁).^51,52 We propagated 600 Wigner-distributed trajectories starting in the S₁ state, with centers at 0.154 449 a.u. for the BLA coordinate and at zero angle for TORS and HOOP, all with zero center momentum. The variances were 0.154 449 a.u. for BLA, 0.183 302 for TORS, and 0.406 143 for HOOP. We used a time-step of 0.01 a.u. a larger time-step gave qualitatively similar results but, for this system, led to a larger energy-conservation violation in CTMQC and also larger jumps in CTMQC-E. Note that Ref. 29 dealt with the energy non-conservation simply by reverting to Ehrenfest dynamics after 120 fs.

Figure 3 shows the quantum yield of the photo-isomerization process, defined as the ratio between the trans and all reaction products, as indicated in the figure. All methods overestimate the quantum yield at short times; Ref. 29 attributed this to the lack of nuclear quantum effects. After about 120 fs, CTMQC-E results oscillate around the exact quantum yield, giving a small improvement over CTMQC, which shows an underestimate, and over EH and SH.

The top panel in Fig. 4 shows the trajectory-averaged total energy as well as the nuclear kinetic energies along the BLA, TORS, and HOOP modes, as compared to EH and SH dynamics and exact calculations. Although not entirely, CTMQC-E strikingly reduces the total energy violation of CTMQC. The residual deviation from constant energy in CTMQC-E is caused by two effects. First, for trajectories for which the denominator in Eq. (17) becomes smaller than a fixed threshold, the accumulated force reverts to its original definition, which does not conserve energy. Second, even when the new definition is used, the quantum momentum computation may revert to the original quantum momentum Q₀ definition when a denominator involved in imposing Eq. (9) becomes too small;^26,42 our redefined force is guaranteed to conserve energy only when used in conjunction with the modified definition Q_m. Considering now the kinetic energy along BLA (the lowest panel), in the exact system, it oscillates around an average that increases in a monotonic way due to the wave packet moving through the conical intersections toward S₀, while the oscillations decrease in amplitude as a result of a loss in vibrational coherence as the wave packet spreads along the TORS degree of freedom. CTMQC captures the initial behavior well but starts increasing at around 80 fs and deviating significantly from the exact calculations; EH and SH, in contrast, over-damp the BLA oscillations and average. On the other hand, CTMQC-E reasonably reproduces the exact kinetic energy along BLA. For the torsional kinetic energy, all MQC schemes yield an overestimate of the kinetic energy from around 40 fs. In CTMQC-E, EH, and SH, this gain in TORS kinetic energy is compensated by a loss in HOOP kinetic energy and a loss in potential energy (not shown here) compared with the exact results. On the other hand, in CTMQC, the kinetic energies of all three degrees of freedom increase simultaneously after 120 fs as a consequence of energy violation. This example appears to be among the most challenging so far for CTMQC for energy conservation, and further work is under way to study this in more detail.

In summary, we have presented a modification of CTMQC, CTMQC-E, that imposes energy conservation over the ensemble of nuclear trajectories through a redefinition of the accumulated force. Much like the modified definition of the quantum momentum that imposes the exact condition of zero population transfer in regions of zero NAC, the new definition involves all trajectories in its computation. It has similarities in form to that used in the independent-trajectory version of CTMQC but results in a more physical condition, being imposed on the ensemble of trajectories rather than on the individual trajectories. Future work includes refinement of the algorithmic implementation to better deal with instabilities arising from small velocities in the denominator of the new definition of the accumulated force and to reduce the jumps occurring in the energy without going to very small time-steps, which may be untenable for larger systems. Generally, for many of the systems where CTMQC has been applied, energy violation is often too small to have much of an effect on its dynamics or physical observables. However as shown here, this is not guaranteed, and future work will also attempt to characterize situations where CTMQC is expected to yield significant energy non-conservation, for which the use of CTMQC-E gives significant improvement.

We thank Dr. Federica Agostini, Dr. David Lauvergnat, and Dr. Lea M. Ibele for their helpful discussions. The financial support from the National Science Foundation under Award No. CHE-2154829; from the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, under Award No. DESC0020044; and from the Computational Chemistry Center: Chemistry in Solution and at Interfaces funded by the U.S. Department of Energy, Office of Science Basic Energy Sciences, under Award No. DE-SC0019394 as part of the Computational Chemical Sciences Program is gratefully acknowledged.