Which Algorithm Best Propagates the Meyer–Miller–Stock–Thoss Mapping Hamiltonian for Non-Adiabatic Dynamics?

A common strategy to simulate mixed quantum-classical dynamics is by propagating classical trajectories with mapping variables, often using the Meyer–Miller–Stock–Thoss (MMST) Hamiltonian or the related spin-mapping approach. When mapping the quantum subsystem, the coupled dynamics reduce to a set of equations of motion to integrate. Several numerical algorithms have been proposed, but a thorough performance comparison appears to be lacking. Here, we compare three time-propagation algorithms for the MMST Hamiltonian: the Momentum Integral (MInt) (J. Chem. Phys., 2018, 148, 102326), the Split-Liouvillian (SL) (Chem. Phys., 2017, 482, 124–134), and the algorithm in J. Chem. Phys., 2012, 136, 084101 that we refer to as the Degenerate Eigenvalue (DE) algorithm due to the approximation required during derivation. We analyze the accuracy of individual trajectories, correlation functions, energy conservation, symplecticity, Liouville’s theorem, and the computational cost. We find that the MInt algorithm is the only rigorously symplectic algorithm. However, comparable accuracy at a lower computational cost can be obtained with the SL algorithm. The approximation implicitly made within the DE algorithm conserves energy poorly, even for small timesteps, and thus leads to slightly different results. These results should guide future mapping-variable simulations.


Introduction
][10] To analyse the complex dynamics observed, accurate and efficient numerical models are of great importance. 3,94][15] This led to the development of approximate classical-like dynamical frameworks, motivated by the classical linear scaling with degrees of freedom (DoF) compared to the exponential scaling of quantum methods. 16Unlike most transformations, discrete quantum DoF do not have an obvious classical counterpart. 16Consequently, approximate methods have been developed to recast the quantum system to look classical whilst retaining some quantum properties, often requiring a compromise between accuracy and cost for large systems.
Many methods exist to incorporate discrete quantum DoF into classical frameworks including Ehrenfest dynamics, 17 approximations to the quantum-classical Liouville equation (QCLE), [18][19][20] the symmetrical quasiclassical windowing method, 21 surface hopping, 22,23 and mapping methods such as methods inspired by the Meyer-Miller-Stock-Thoss (MMST) mapping, [24][25][26][27][28] and spinmapping, 29,30 including the mapping approach to surface hopping (MASH). 31The simplest nonadiabatic method, proposed by Mott (1931), evaluates the quantum electronic dynamics along the classical path of the nuclei, known as the 'classical path' approach. 16,18,19,32However, this does not include the 'back reaction' on the nuclei. 16urface hopping models, originally developed by Tully  et al. (1971), propagate along one adiabatic surface before 'hopping' to another. 16,22,23,33,34'Hopping' is possible at any point along a trajectory, not just where surfaces cross, destroying the state coherence. 23,26 this article, we focus on mapping approaches, utilised in various semiclassical methods, where the discrete quantum DoF are mapped onto continuous classical DoF propagated by classical mechanics.The Meyer-Miller mapping developed in 1979, 25 and later put on a rigorous footing by Stock and Thoss in 1997, 26 constructs a set of classical variables for the discrete electronic DoF and propagates them with the nuclear DoF using an effective Hamiltonian, the MMST Hamiltonian. 25,26,35he electronic dynamics are consistent with the timedependent Schrödinger equation, and the force exerted on the nuclei is given by the instantaneous values of the electronic variables. 25,35The mapping introduces electronic position and momenta, sometimes expressed as action-angle variables, to describe the nuclear motion on coupled potential energy surfaces. 16,26This approach maps the time-dependent Schrödinger equation for an Nlevel system to a classical analog of N coupled harmonic oscillators following Hamilton's equations of motion. 16,26he recently derived spin-mapping approach employs a different mapping formalism but gives a Hamiltonian almost identical to that of the MMST methods. 29,30,36The MMST Hamiltonian algorithms compared in this work are thus also applicable to spin-mapping methods.Compared to MMST-based methods, spin-mapping uses a different value of the so-called zero-point energy (ZPE) parameter, introduced by Stock and Müller as a fitting parameter to mitigate ZPE-leakage. 37,38Spin-mapping leads to a ZPE as a function of the number of states, with values close to what was previously found optimal when tuned as a free parameter.A partially linearized spin-mapping method has been found to improve accuracy compared to the fully linearized mapping, 39,40 in particular for spectroscopy. 41ne application of the MMST Hamiltonian is to calculate dynamical properties at thermal equilibrium by approximating equilibrium time-correlation functions where the partition function Z is given as and the quantum Boltzmann operator at inverse temperature β = 1/k B T is e −β Ĥ .Semiclassical methods can be used to calculate quantum time-correlation functions, often by utilising an 'Initial-Value Representation' (IVR), resulting in a phase-space integral. 1,26,42The semiclassical phase-factor makes correlation function convergence challenging. 1,27,438][49][50][51] However, semiclassical methods often fail at describing nuclear quantum effects. 525][66] Specifically, RPMD provides, for a single (adiabatic) potential, an approximation to Kubo-transformed time correlation functions while preserving the Boltzmann distribution and is exact in the shorttime and classical limits. 58,59,67The Kubo-transformed correlation functions allow short-time quantum effects to be included, although the long-time quantum coherence effects are neglected.4][75][76][77] However, none of these methods alone fulfils the three important limits of; replicating Rabi oscillations in the uncoupled case, preserving the quantum Boltzmann distribution, and reducing to classical dynamics in the adiabatic limit. 13,29,78hile a recently developed ellipsoid spin-mapping fulfils all these limits, its mean-field dynamics was found to be often less accurate for short times than the original spinmapping. 78Further work is still required to find an accurate trajectory-based approach to replicate the true quantum dynamics.
Although mapping Hamiltonians have become increasingly popular for simulating non-adiabatic dynamics, numerical integration of the equations of motion is not straightforward due to coupling between the electronic momenta and nuclear positions.Over the years various algorithms have been suggested to solve this problem, [1][2][3] but as far as we are aware there has been no published computational comparison of MMST algorithms to determine their properties and accuracy.2][3] Other algorithms exist for this problem, such as Runge-Kutta or the Adams-Bashforth Predictor-Corrector algorithm which are known to be non-symplectic.We note that one cannot directly use velocity-Verlet due to coupling between the nuclear positions and the mapping momenta.To avoid integrating stiff equations of motion, Wang et al. suggested an intelligent canonical transformation. 27However, even with this transformation, these approaches are not ideal for propagating mapping variables and will still require short time steps.We will not consider these algorithms further and instead focus on algorithms which attempt to propagate the mapping variables exactly for an arbitrary timestep.
We seek to compare the algorithms for the simplest possible system for which they can all be compared on an equal footing, and for which there already exists results in the literature for qualitative comparison.To this end, we compute position and state autocorrelation functions for a two-state linear vibronic potential with the MMST Hamiltonian (corresponding to the single-bead limit of NRPMD) using these three algorithms. 53As many of the methods discussed share similar Hamiltonian forms to the MMST, including spin-mapping, Ehrenfest and some surface hopping models, the results are widely applicable and should inform future computational studies using mapping variable methods.
The article is structured as follows.In section 2 we provide theoretical background for the three algorithms investigated.In section 3 we investigate the symplecticity, satisfaction of Liouville's theorem, accuracy of trajectories and correlation functions, computational cost and energy conservation.We conclude in section 4.

Background Theory
Here, we present the algebraic forms of the three algorithms on an equal footing, such that we can compare properties including the accuracy, symplecticity and if they satisfy Liouville's theorem.
Note that a brief theoretical analysis of the MInt and SL algorithms, determining the symplecticity in particular, was previously carried out in Ref. [1].We extend the work in Ref. [1] to a more thorough computational investigation where we compare the algorithmic performance utilising a model for which there exists results in the literature, and include the DE algorithm.As far as we are aware, the symplecticity has not been rigorously determined for the DE algorithm.

Symplectic Integrators
The Hamiltonian for an N-level electronic system in the diabatic representation is where the nuclear position and momenta are x and p respectively, V(x) is an N × N diabatic electronic potential energy matrix in the basis of the electronic states, |n , and µ µ µ is a diagonal matrix of nuclear masses.The classical MMST mapping Hamiltonian in the diabatic representation is 25,26 (4) where X and P are the electronic position and momenta respectively.Classical evolution under this Hamiltonian refers to the equations of motion 16,26 that constitute a time-dependent Hamiltonian system where z = [x T , X T , p T , P T ] T , the ∇ z operator contains the partial derivatives, J is the structure matrix and I/O are the identity and zero matrices respectively. 1,25,79A Hamiltonian integrator is said to be symplectic if it fulfils the condition 1,79,80 where M is the monodromy matrix.The monodromy matrix is a matrix of differentials that expresses how the time-evolved phase-space variables depend on the initial phase-space variables 81 The monodromy matrix is computed for each timestep and multiplied with the previous timestep matrices to obtain the monodromy matrix for the overall trajectory If the Hamiltonian is split, H = H 1 + H 2 , then the monodromy matrix is calculated for each timestep according to the algorithmic flow map. 79Calculation of the monodromy matrix is not needed for algorithm propagation, only to determine symplecticity and if the theoretical framework requires it.For example, in some SC-IVR methods the phase-factor can be calculated using elements of the monodromy matrix, where utilising a symplectic algorithm can improve the stability, partially mitigating the notorious 'sign' problem. 1,82,83e symplecticity criterion, Eqn.(8), is a much stricter condition than the conservation of volume phase-space (Liouville's theorem), which only requires the monodromy matrix determinant to be unity. 1,80,84Volume phase-space preservation is a consequence of symplecticity, but this relationship does not necessarily hold the other way around, 79 i.e. volume phase-space preservation is necessary but not sufficient for symplecticity.Instead, symplectic integrators can arise from exact time-propagation of a Hamiltonian system. 79Splitting the Hamiltonian into subevolutions will also result in a symplectic integrator, provided that each sub-evolution is the exact timepropagation of the relevant sub-Hamiltonian. 1,79For example, velocity-Verlet is a classical symplectic algorithm as splitting results in two sub-Hamiltonians that are independent of each other, such that both can be integrated exactly. 80The MMST Hamiltonian contains coupling between the electronic momenta and nuclear position, making symplectic time-evolution of the equations of motion challenging and in general, results in nonlinear dynamics. 1,26Symplectic integrators are an advantage for Hamiltonian integration as they have little to no energy drift with time and tend to be more stable at long simulation times. 1 We also note that the Cayley transform can improve the symplectic stability of algorithms with no additional algorithmic complexity, computational cost and can be implemented for any path-integral based scheme. 85

Algorithmic Overview
The MInt algorithm, by Church et al. (2018), was developed to help extend the MQC-IVR method to simulate nonadiabatic dynamics using the MMST mapping. 1 The name arose as the MInt algorithm exactly solves the Momentum Integral with time and as we will show, is the only known symplectic algorithm propagating the MMST Hamiltonian. 1 The MInt algorithm splits the MMST Hamiltonian into two sub-Hamiltonians, each of which is propagated exactly Exact evolution of the sub-Hamiltonians results in symplecticity, and the sub-evolution of H 1 is split into two half timesteps to improve the time-order error, such that the algorithm is at least a second-order method. 1,79Hamilton's equations of motion are obtained for H 1 and H 2 , with the latter being more complicated due to coupling between the nuclear and quantum DoF. 1 The MInt algorithm can be used on any Hamiltonian containing a sum of Meyer-Miller-like terms and has algebraically been shown to be symplectic, symmetric, second-order in time and timereversible. 1Fortran code of the MInt algorithm is available in the SC-IVR package on the Ananth Group website. 86ecently, the MInt algorithm was utilised by Gardner et al. in NQCDynamics.jl,a Julia package for condensed phase nonadiabatic quantum dynamics. 87To be able to compare the three algorithms on an equal theoretical footing, we use the D5 form of the MInt algorithm where instead H 2 is split into two half-timesteps 1 where Ψ MInt H,∆t is the approximate flow map comprised of exact evolutions of the relevant sub-Hamiltonian, Φ H,∆t , and propagation is done from right to left.In our notation, Φ refers to exact evolution and Ψ refers to approximate evolution which may or may not be comprised of exact subevolutions, consistent with the notation of Leimkulher and Reich. 79e Split-Liouvillian (SL) algorithm, by Richardson et  al. (2017), uses the Liouvillian formalism to construct an integrator. 3The Liouvillian operator can be generated from the Hamiltonian using a Poisson bracket where the solutions to time-evolution are of exponential form. 84We follow the convention in the Matsubara dynamics article 59 and by Zwanzig, 88 and define the Liouvillian to be real with no factor i (the imaginary unit).The MInt algorithm D5 form, Eqn.(12), is equivalent in where evolution requires propagation of the Liouvillians from the right to the left and L 1/2 are the Liouvillians of H 1/2 respectively.The SL algorithm further symmetrically splits L 2 into an electronic (el) and nuclear (p) contribution.Hence, the flow map for this algorithm is where are the electronic and nuclear Liouvillians of H 2 respectively, where i is the nuclear index and j is electronic index. 1,3This uses the approximation e L 2 ∆t = e (L p +L el )∆t e L p ∆t e L el ∆t , which is valid in the limit ∆t → 0. 3 Hence, the flow map can also be represented as where Ψ SL a/b represents the approximate propagation of H 2 .Due to symmetric splitting the two half timesteps are not equal and are labelled with the superscripts a for e L p e L el and b for e L el e L p .For an integrator to be symplectic, it is sufficient but not necessary that the integrator is a sequence of exact sub-Hamiltonian evolutions.Evolution under Liouvillians is not guaranteed to be symplectic but will be if each Liouvillian corresponds to exact propagation of a Hamiltonian or sub-Hamiltonian. 79For an arbitrary timestep, the splitting of H 2 into an electronic and nuclear contribution results in holding the electronic position and momenta constant while propagating the nuclear momenta.Hence, the propagation of H 2 is no longer exact or guaranteed to be symplectic. 1 Church et al. (2018) derived the SL monodromy matrix and using the symplecticity criterion, Eqn. ( 8), confirmed that the SL algorithm is not symplectic. 1 They stated that there is likely to be an energy drift associated but this may be small if the adiabatic states are close in energy and weakly coupled. 1 However, the algebraic proof does not indicate if the SL algorithm will result in unreasonable dynamics such as a very large energy drift when compared to the MInt for a given timestep.Kelly et al. (2012) have also developed an integrator that we refer to as the DE algorithm due to the approximation needed to obtain the final algebraic form.The DE algorithm also splits the Hamiltonian into H 1 and H 2 , where the motivation is to enhance stability and minimise the difference between the exact and approximate dynamics. 2The algorithm uses MInt-like equations but makes an implicit approximation equivalent to assuming that the diabatic potential matrix eigenvalues are degenerate, 89 simplifying the calculation of the potential matrix derivative, making the DE algorithm unlikely to be symplectic. 1,2The flow map cannot easily be expressed in a Liouvillian form but can be written as where The degenerate eigenvalue approximation has previously been utilised in PBME simulations of coherent dynamics in photosynthetic systems. 90However, this is known to cause the state populations to deviate from the exact results for systems with large energy biases. 892][3] The DE algorithm has a similar form to the SL algorithm, but the DE algorithm does not hold the electronic position and momenta stationary while propagating the nuclear momenta and has a different diabatic potential matrix differential with respect to nuclear position. 2,3As far as we are aware, the DE algorithm monodromy matrix has not been algebraically determined, tested for symplecticity or Liouville's theorem.

Propagation of H 1
2][3] All three algorithms propagate H 1 in the same way using Hamilton's equations of motion [1][2][3] where m is the nuclear mass, such that integration provides the propagation equations We define our monodromy matrix as such that the monodromy matrix for the propagation of H 1 is a triangular matrix of the form where 0 = [0, 0] T and the determinant is unity, The propagation of H 1 is symplectic for all three algorithms as M T H

Propagation of H 2
Here, we present and compare the propagation of H 2 for the three algorithms.The diabatic potential matrix is split into a state-independent term, U, and a traceless statedependent matrix, Ṽ, such that V(x) = U(x) + Ṽ(x).This is as the DE algorithm requires a traceless state-dependent matrix, previously being seen as an advantage as it renders the dynamics invariant to any constant shift of the coupling potential. 2,91This gives where X = [X 1 , . . .][3] The propagation of H 2 is split into two half-timesteps that sandwich Φ H 1 , allowing fair algorithmic comparison with the original SL form and the equivalent MInt form.We have swapped the order of H 1 and H 2 to test the DE algorithm and in Figs. 5 and 6 in the SI, we show that this has no effect on the symplecticity or energy conservation.

The Momentum Integral (MInt) Algorithm
The MInt algorithm exactly propagates H 2 with time, taking into account the electronic dependence of the nuclear momentum. 1Using Hamilton's equations of motion for half a timestep, following the approach of Church et al. we derive propagation equations where the prime denotes the derivative with respect to x and the matrix dependence on x has been dropped for simplicity.To find solutions through integration, the MInt algorithm relies on the fact that Ẋ and Ṗ are independent of p but ṗ is dependant on X and P. 1 Therefore, solving for X(t + ∆t 2 ) and P(t + ∆t 2 ) and substituting into Eqn.(25b) allows p(t + ∆t 2 ) to be found.The electronic position and momenta are found through integration resulting in which can be recast to be entirely real by diagonalizing Ṽ into eigenvectors, S, and a diagonal eigenvalue matrix, Λ Λ Λ, such that S T ṼS = Λ Λ Λ. 1 This results in for the propagation of X and P, where C and D are, 1 By defining the derivative of the potential in the adiabatic basis to be G = S T Ṽ S and inserting SS T = I identities into the integration of Eqn.(25b), an intermediate equation for the propagation of p can be found which can be solved by element-wise integration of Defining where λ nm = (Λ Λ Λ) mm − (Λ Λ Λ) nn , such that where E is symmetric and F is skew-symmetric.The nuclear propagation, where x is unchanged, is therefore The extended form for multiple states can be found in Appendix B of Ref. [1].The MInt monodromy matrix for the propagation of H 2 satisfies Liouville's theorem and can be found by defining such that the monodromy matrix is 1 The electronic propagation resulting from L el is equivalent to the MInt algorithm, Eqn.(27), and leaves the nuclear variables unchanged. 1L p results in the following propagation of p and leaves all other variables unchanged The difference between the two algorithms arises in the nuclear momentum, where due to the SL symmetric propagation of H 2 , the electronic variables are always the half-timestep evolved values for propagation of p. 1,3 The SL monodromy matrices for the propagation of H 2 obey Liouville's theorem and are where Neither M p and M el satisfy the symplecticity criterion individually, seen in Appendix A. 1 The combination of M p and M el was shown by Church et al. to only be symplectic in the ∆t → 0 limit, not for an arbitrary timestep.

The Degenerate Eigenvalue (DE) Algorithm
The DE algorithm defines Ṽ to be traceless, such that the last term of H 2 is ignored. 2 Here, we will re-frame the algorithm into a MInt-like form using vector notation and outline the degenerate eigenvalue approximation used.We will also investigate the symplecticity and satisfaction of Liouville's theorem through defining the monodromy matrix, which as far as we are aware has not been done before.
The electronic Hamiltonian equations are the same as above for the MInt and SL algorithms and solved to provide the same overall propagation equations, Eqn.(27).The nuclear momentum propagation equation is To solve Eqn.(39), we define an integral, A, that requires the degenerate eigenvalue assumption to arrive at the final form by Kelly et al. (2012) This can be written in a MInt-like form such that using Eqn.( 26) A transformation is defined into the adiabatic basis using the following overlined variables 2 X = S T X , P = S T P, such that, conversion into the adiabatic basis utilising the decomposition of Ṽ results in where Ṽ is determined by We can arrive at the form shown in Ref. [ 2] by making the approximation that the eigenvalues are equal, Λ Λ Λ εI.Therefore, the differential of the eigenvalues can be approximated as Λ Λ Λ ε I such that the derivative of Ṽ in the adiabatic basis is obtained through differentiating the identity, I = (SS T ) = S T S + S T S = 0, Inserting this into Eqn.(45) results in Hence, the integral below is obtained The final form seen in Ref. [2] has an additional − ∆t 4 Tr[Λ Λ Λ ] term when written in matrix form, however, this is just a dummy term.The fact that Ṽ is traceless requires that the sum of the eigenvalues is zero, therefore, Λ Λ Λ and Λ Λ Λ are traceless.The overall nuclear propagation in the diabatic basis is then To obtain the final form by Kelly et al. (2012), we make the DE approximation to simplify Eqn.(47). 2,892][3] It can be seen that in the case where the DE approximation holds such that, Ṽ = SΛ Λ Λ S T , the propagation of p would be similar to the SL algorithm.However, the DE algorithm uses the inital electronic variables when the SL algorithm uses the half-timestep evolved values.Due to the similarities prior to the DE approximation being made, the DE algorithm is unlikely to be symplectic.To rigorously check the symplecticity, we derive the monodromy matrix by defining such that The determinant of M DE can easily be shown to be unity, satisfying Liouville's theorem.However, the approximate propagation of H 2 is not symplectic as M DE does not satisfy the symplecticity criterion, shown in Appendix A. The MInt algorithm is the only algorithm considered here that is both symplectic and satisfies Liouville's theorem. 1 Below, we present a table of the theoretical results for easy algorithmic comparison.1: Summary of the theoretical results obtained here and Ref. [ 1].The H 2 electronic propagation for all algorithms H 2 is equivalent.The approximations made in the SL and DE algorithms result in inexact nuclear momentum propagation.However, the SL algorithm is exact in the ∆t → 0 limit.

Results and Discussion
The algorithms discussed here have been utilised in the literature but not for the same system to allow direct comparison.We seek to test the algorithms computationally on an equal footing using the same system.Hence, we use a simple two-state linear vibronic potential, also known as a double well potential, discussed in this section with the MMST Hamiltonian (corresponding to the single-bead limit of NRPMD) to test the symplecticity, satisfaction of Liouville's theorem, energy conservation and accuracy of correlation functions.The models defined here also allow qualitative comparison with literature.

Theoretical Models
To compare the approximate dynamics produced by the MInt, SL and DE algorithms, we use the three potential models introduced in Ref. [ 53].For these models the potential diabatic matrix is (52)   such that splitting to obtain a traceless Ṽ gives Reduced units are used where m = h = ω = 1, so energy is measured in units of the frequency, ω.The models represent bound potentials, defined in Table 2, where ∆ is the electronic coupling, 2α is the energy bias between the potential energy surfaces (the asymmetry) and κ is the vibronic coupling, chosen to be Models 1-3 were utilised to compute correlation functions using the following distribution, ρ, in the N-bead form in Refs.[ 53, 3].Here, for simplicity, we use the single-bead form such that the partition function is where W = P T MXX T MP and M = e −β Ṽ/2 .Note that W is positive definite for the single bead case.Monte Carlo importance sampling is utilised as detailed in Appendix B. Time-independent equilibrium properties are calculated, for example, the population of state n 3,53 The derivation is similar to that of Eqn.(B.4) in Appendix B by summing over all indices except n, 3,53 resulting in where, This approach for obtaining population information is only valid at t = 0, whereas the electronic populations can be found at time t using 3 Correlation functions are calculated through finding an approximation to Eqn. (1), where for the position autocorrelation function Â = B = x and for the population autocorrelation function Â = A n and B = B n , where W (X i (0), P i (0)) and initial conditions are sampled from Eqn. (B.14) for J trajectories.The correlation functions are then averaged over J trajectories where index i refers to the ith trajectory.We can directly compare the three algorithms and qualitatively compare our results with Ref. [ 53].We sample the same distribution, with our method outlined in Appendix B, as in Ref. [ 53].However, we choose to split the potential matrix such that Ṽ is traceless as the DE algorithm requires this.In Ref. [53], the matrix was instead split such that the lowest eigenvalue of Ṽ is zero, which in general results in a non-zero trace.
To ensure we can compare all three algorithms for the same Hamiltonian, we have used the traceless form leading to small quantitative differences between the results here and in Ref. [53].

Algorithmic Properties
First, we consider the symplecticity and conservation of Liouville's theorem using Model 1.
Church et al. determined that the MInt algorithm is symplectic whilst the SL algorithm is not. 1 In Appendix A, we algebraically show that the DE algorithm is not symplectic.
To numerically determine the symplecticity, we define an error matrix, E r , to be where for a symplectic integrator the elements of E r , a i j , will all be zero. 92The Frobenius Norm is used to track the size of E r where the matrix size is n × n. 92 To average over many trajectories, we weight by where i refers to the trajectory index.To determine if Liouville's theorem is satisfied, we evaluate which will be zero if it is satisfied. 1,80By squaring the deviation, we ensure no error cancellation when averaging over trajectories.In Fig. 1a, the logarithmic plot of ||E r || F against time for Model 1 with ∆t = 0.1 can be seen.The MInt algorithm (cyan) remains below 10 −12 for the entire simulation time and is therefore symplectic, with a small build up of floating-point errors that arise in numerical calculations, in agreement with the literature. 1 The SL algorithm (purple) increases rapidly to around 10 −2 and continues to increase indicating that it is not symplectic.The DE algorithm (red) is the least symplectic, being on the order of 1 by the end of the simulation time.Although the SL and DE algorithms are not symplectic, all three algorithms satisfy Liouville's theorem and conserve volume phase-space, as seen in Fig. 1b.We believe that the very slight increase seen for the MInt algorithm is due to additional floatingpoint error accumulation arising from the more complicated propagation equations.We now look at energy conservation, Fig. 2, where (a) depicts the energy of a single trajectory and (b) averages the energy conservation criterion, over trajectories until convergence was observed.We calculate the energy, ε, by evaluating the MMST Hamiltonian, Eqn. 4, at each timestep.Under perfect energy conservation, the criterion should be zero.For a single trajectory, we observe that the DE algorithm has much larger oscillations and does not conserve energy well when compared to the MInt and SL algorithms.In Fig. 2b, we average over many trajectories and consider different timestep sizes.It is seen that the MInt and SL algorithms are second order, as changing the timestep by a factor of 10 increases the criterion by ∼ 10 4 as expected. 79his agrees with the algebraic determination of the order by Church et al.. 1 However, the DE algorithm stays the same magnitude for ∆t = 0.01 and ∆t = 0.1, appearing to not be affected by the smaller timestep.The very coarse ∆t = 1.0 is so large that it breaks the trends.The DE algorithm's poor energy conservation is due to the discarded terms when the DE approximation is made; when these terms are large, the propagation of the nuclear momentum is affected significantly which then affects all other variables through the propagation equations.The MInt and SL energy conservation is very similar although the SL algorithm is seen to have the smallest energy fluctuations throughout and appears to have a longer period of oscillation for Model 1.This is surprising as one would expect a symplectic algorithm to have better energy conservation compared to a non-symplectic algorithm.For Models 2 and 3, we observe that the MInt and SL algorithms have almost identical energy conservation whilst the DE algorithm is worse, seen in Figs. 1 and 2 in the SI.

Correlation Functions
The nuclear position and electronic population autocorrelation functions, Eqn.(61), were calculated for the three models.The fast oscillations of C11 in Fig. 3a indicate that the strong electronic coupling is close to the adiabatic limit and the dynamics tend towards standard RPMD as U is approximately harmonic. 53For Models 2 and 3, the weaker electronic coupling reduces the electronic oscillation frequency, producing curves that deviate from the adiabatic result. 53In Model 2, the equilibrium population of the first electronic state is quickly lost to the lower energy second state, indicating that the system is almost always on one diabatic surface.The correlation functions obtained using the MInt and SL algorithms qualitatively replicate the dynamics expected from the single-bead calculation in Fig. 1 of Ref. [ 53] for all models.Comparing the three algorithms tested in Fig. 3 provides the interesting discovery that the DE algorithm is very accurate for Model 1, despite the lack of energy conservation.This is likely due to averaging with fast electronic oscillations providing the correct convergence.However, for the other models the DE algorithm predicts the same initial drop in the electronic population autocorrelation functions as the MInt and SL algorithms, but starts to deviate from the expected behaviour after the minima.The correlation functions have the same shape indicating a systematic error that likely arises due to the DE approximation.This assumes the offdiagonal elements of Ṽ in the adiabatic basis, G, are zero and that the diagonal elements are equal, which is not the case for the SL and MInt algorithms. 1 Whilst testing energy conservation, we observed that the DE algorithm has more trajectories with poor energy convergence for models with weaker electronic coupling and was particularly poor in the intermediate regime given by Model 3. The MInt and SL algorithms produce the same correlation functions for all the models tested with small timesteps.When testing different timesteps, as seen in Appendix A, both the MInt and SL algorithms are tolerant of a coarse timestep.For ∆t = 1.0, we observe that the MInt and SL correlation functions start to differ, with the MInt being closer to the small timestep results.However, the electronic oscillations are not captured well due to aliasing.The MInt and SL algorithms are limited by the model used rather than the algorithmic accuracy.

Computational Results
MInt SL DE Satisfies Liouville's Theorem Symplectic Energy Conservation Good Good Poor Correlation Function Accuracy Good Good Poor Table 3: Summary of the computational results, which are consistent with the theoretical results in Table 1.
In Table 3, we provide an overview of the properties tested here.3]53 The lack of symplecticity may result in an upwards energy drift at long simulation times, although this was not observed within the short simulation time tested here.
The SL algorithm approximates that the electronic variables can be held still while the nuclear momentum is propagated. 3The DE algorithm assumes that the eigenvalues of the potential matrix are equal, such that the potential derivative in the adiabatic basis is approximated as S T Ṽ S = Λ Λ Λ . 2 This results in inaccurate nuclear trajectories that leads to poor energy conservation.However, all three algorithms obey Liouville's Theorem, preserving volume phase-space throughout the trajectories.The MInt and SL have similar energy conservation and are both second-order algorithms, where in some cases the SL has better energy conservation. 1,79he correlation functions produced for the MInt and SL algorithms are very similar to those computed in Ref. [53], with the small differences arising from the different choice of splitting the potential.The DE algorithm converges to a different result for models with weaker coupling, being a good approximation only near the adiabatic limit.
Table 1 in the SI presents timings of the computational algorithms.One should note that in this case, the bottleneck is the evaluation of the algorithm itself rather than the calculation of the potential matrix.When calculating the monodromy matrix, we observed that the SL algorithm has the lowest computational cost followed closely by the MInt, with the DE algorithm taking the longest time to run.Without calculating the monodromy matrix, the SL is significantly faster and the MInt and DE take similar times to run.

Conclusions
In this article we have tested symplecticity, Liouville's theorem, energy conservation and computed correlation functions using the MInt, SL and DE algorithms for a range of model parameters and timesteps.We find that the computational results agree with our theoretical predictions.If symplecticity is required, for accurate MMST Hamiltonian dynamics with little energy drift, the MInt algorithm should be used.As far as we are aware, the MInt is the only known symplectic algorithm for a general form of the MMST Hamiltonian.However, even though the SL algorithm is not formally symplectic with a finite timestep, it becomes exact in the limit of an infinitesimal timestep.In our tests, it gave comparable accuracy to the MInt algorithm, but at a lower computational cost.We would not recommend the DE algorithm for these models, as it breaks energy conservation and introduces errors into the results.This indicates that for the models used here, the degenerate eigenvalue approximation is not valid.
Further work includes integrating the Cayley transform to extend these findings to stable NRPMD simulations. 85dditionally, one can apply the MInt algorithm to related dynamical methods such as forward-backward (FB)-IVR. 51,82

A Symplecticity
The symplecticity criterion can be algebraically determined for the three algorithms, where the MInt and SL have already been derived by Church et al. (2018).It is sufficient to say that under the splitting of the Hamiltonian into H 1 and H 2 , the symplecticity criterion becomes 1,79 For example, following the method by Church et al., M H 1 is Eqn.( 23) for all three algorithms.This can be derived from the propagation equations, Eqn. ( 21), where Evaluating the symplecticity criterion results in such that the propagation under H 1 is symplectic. 1 The symplecticity for H 2 is derived separately for all three algorithms due to the different nuclear propagation equations, where the MInt and SL algorithms had previously been derived by Church et al..For the MInt algorithm, M H 2 is Eqn.(35) where the symplecticity criterion becomes where Church et al. previously determined that h ≡ 0 and j ≡ 0 ∀ X, P.
For the SL algorithm, propagation of H 2 is split further into a nuclear and electronic contribution for which the monodromy matrices are Eqn.(38).Evaluating the symplecticity criterion for M p M where Eqn.(A.5) has been used in conjunction with the fact that h ≡ 0 and j ≡ 0 ∀ X, P. 1 Church et al. note that a, e = 0 T so evolution under L el and L p is not symplectic.The combined evolution under L el and L p does not provide error cancellation that restores symplecticity.We find that numerically the difference between M H 2 and M el M p is on O(∆t) 2 , in agreement with the literature. 1 The combination of M p M el also provides the same result.This means the propagation of M p M el will be symplectic in the ∆t → 0 limit, but will not be for an arbitrary timestep.
For the DE algorithm, the monodoromy matrix, M DE , for the approximate propagation of H 2 has been derived for the first time as Eqn.(51).Using Eqn.(34) with h ≡ 0 and j ≡ 0 ∀ X, P, the symplecticity criterion becomes which will not be symplectic for the same reasoning as for the SL algorithm.
To derive the order of the difference between M H  The Frobenius Norm of the error matrix averaged until convergence with the SL (purple), the MInt (cyan), and the DE (red) algorithms using timesteps of ∆t = 0.1 (solid), ∆t = 0.01 (dotted) and ∆t = 1.0 (dashed).The MInt algorithm is seen to be symplectic and does not change with timestep.The SL algorithm is seen to be second order with respect to time and the DE algorithm is zero order, in agreement with the theoretical predictions.
We find that, numerically, when considering the whole propagation of H using the SL algorithm, the initial difference in symplecticity between the MInt and SL algorithms scale on order of O(∆t) 3 .This is due to some error cancellation with the symmetric splitting of H 2 .Hence, when propagating for a given length of time, the order is O(∆t) 3

B Sampling
Following a similar derivation of the partition function as in Refs.[ 3] and [ 53], the sampling method used in this work is derived.Instead of formulating in the limit of an infinite number of ring-polymer beads, we derive the partition function in the single-bead case.The wavefunctions of the singly-excited oscillator (SEO) states are known in the position and momentum bases 3 This means that W = (X T MP) 2 , which will always be positive.Sampling from Gaussian distributions for electronic variables where µ = 0 and σ = 1/ √ 2, and classical distributions for nuclear variables where µ x = 0 and σ x = 1/(mω 2 β ), and µ p = 0 and σ p = m/β such that the sampled distribution is Hence, we need to correct observables calculated to average over the required distribution.We can thus define the average of an observable over the distribution, Â, that may be a function of initial position and momenta i.e, Â(X i , P i , x i , p i ) as × ρ samp dxdpdXdP. (B.16) The integral over the sampled distribution can be replaced with a sum over the number of samples taken, J, The partition function can be calculated as

Figure 1 :
Figure 1: (a) The Frobenius Norm of the symplecticity error matrix and (b) the determinant criterion as a function of time using Model 1 and ∆t = 0.1, averaged over a million trajectories using the SL (purple), the MInt (cyan) and the DE (red) algorithms.In (a) the MInt algorithm is seen to be symplectic whereas the DE and SL are not, whereas in (b) all algorithms satisfy Liouville's theorem.

Figure 2 :
Figure 2: The energy conservation for Model 1 with (a) a single trajectory and ∆t = 0.1 (solid) and (b) averaged using ∆t = 0.1 (solid), ∆t = 0.01 (dotted) and ∆t = 1.0 (dashed) with the SL (purple), the MInt (cyan), and the DE (red) algorithms.The DE algorithm has the worst energy conservation and, for this system, the MInt has slightly worse conservation than the SL.The SL and MInt algorithms show improved energy conservation upon decreasing the time step, unlike the DE.

Figure 3 :
Figure 3: The nuclear position, Cxx (t), and electronic population, C11 (t), autocorrelation functions for (a) Model 1, (b) Model 2 and (c) Model 3 using ∆t = 0.1 with the DE (red solid), SL (purple dashed) and MInt (cyan dotted) algorithms.The SL and MInt give identical results, but the DE algorithm deviates from them in (b) and (c).

Figure A. 1 :
Figure A.1:The Frobenius Norm of the error matrix averaged until convergence with the SL (purple), the MInt (cyan), and the DE (red) algorithms using timesteps of ∆t = 0.1 (solid), ∆t = 0.01 (dotted) and ∆t = 1.0 (dashed).The MInt algorithm is seen to be symplectic and does not change with timestep.The SL algorithm is seen to be second order with respect to time and the DE algorithm is zero order, in agreement with the theoretical predictions.
/∆t = O(∆t) 2 , seen in Fig. A.1.For the DE algorithm, the initial difference in symplecticity with the MInt algorithm scales on O(∆t).Therefore, for a given length of time, the order is O(∆t)/∆t = 1.Fig. A.1 corroborates this as the DE algorithm appears to be on zero order for all the timesteps tested.The symplecticity does not improve with the smaller timesteps.

Table 2 :
Model 2 has a strong energy bias and the system is in the inverted Marcus regime.Model 3 is a challenging intermediate regime in which the timescales of the electronic and nuclear dynamics are similar.Values for the potential matrix constants for Models 1-3 from Ref.[53].
1. Model 1 represents strong electronic coupling, where the nuclear dynamics occur on a longer timescale than the electronic oscillations and nuclear motion occurs in a mean field of the diabatic surfaces.