Recovering Marcus Theory Rates and Beyond without the Need for Decoherence Corrections: The Mapping Approach to Surface Hopping

It is well-known that fewest-switches surface hopping (FSSH) fails to correctly capture the quadratic scaling of rate constants with diabatic coupling in the weak-coupling limit, as expected from Fermi’s golden rule and Marcus theory. To address this deficiency, the most widely used approach is to introduce a “decoherence correction”, which removes the inconsistency between the wave function coefficients and the active state. Here we investigate the behavior of a new nonadiabatic trajectory method, called the mapping approach to surface hopping (MASH), on systems that exhibit an incoherent rate behavior. Unlike FSSH, MASH hops between active surfaces deterministically and can never have an inconsistency between the wave function coefficients and the active state. We show that MASH not only can describe rates for intermediate and strong diabatic coupling but also can accurately reproduce the results of Marcus theory in the golden-rule limit, without the need for a decoherence correction. MASH is therefore a significant improvement over FSSH in the simulation of nonadiabatic reactions.

Introduction: Under the Born-Oppenheimer approximation, one assumes that electronic motion is fast compared to nuclear motion and is therefore adiabatically separated.The resulting picture of nuclei moving on a single adiabatic potential energy surface forms the basis of our modern understanding of molecular structure and dynamics.Despite its great success, there are many important molecular processes for which the Born-Oppenheimer approximation is not valid.6][17][18][19] Fortunately, however, the relatively high mass of atomic nuclei means that it is often a reasonable approximation to treat them as classical particles with well-defined positions and momenta.In 1990 Tully proposed what has become the most widely used of such 'mixed quantum-classical' methods for simulating nonadiabatic processes, known as fewest-switches surface hopping (FSSH). 20Within FSSH, the nuclei predominantly move under the force of a single adiabatic potential energy surface, with occasional stochastic hops between the surfaces.The probabilities for these hopping events are determined based on the evolution of the electronic wavefunction under the time-dependent Hamiltonian generated by the nuclear trajectory.
2][23][24] The result is a deviation between the number of trajectories on each surface and the wavefunction coefficients, which can therefore be referred to as an inconsistency error.][25][26][27][28][29] Due to their ad-hoc nature, decoherence corrections are not guaranteed to consistently improve the results of a calculation. 30However, one area where they have been shown to be essential is for processes, such as electron transfer, which involve slow population transfer in strongly nonadiabatic systems (weak diabatic coupling, ∆).2][33][34][35] This was explained in terms of repeated crossings of the nonadiabatic coupling region, leading to a build up of the inconsistency error. 31ecently, an alternative to FSSH has been derived known as the mapping approach to surface hopping (MASH). 36ASH was designed to offer the best of both worlds between surface hopping and mapping approaches, such as the Meyer-Miller-Stock-Thoss mapping 37,38 and spin mapping. 39,402][43][44][45][46] Tests against exact results for the Tully models, a series of spin-boson models, as well as 3-mode and 24-mode vibronic models of pyrazine have shown that the results of MASH are generally as good or better than FSSH for an equivalent computational cost. 36erhaps most interesting are the results for the spin-boson model where the system crosses the coupling region many times during the dynamics.One might have expected that decoherence corrections were necessary to improve upon the FSSH results.However, MASH shows a significant improvement even without the addition of decoherence corrections.This raises the question, how well will MASH perform in systems exhibiting slow population transfer with weak diabatic coupling where the errors of FSSH are known to be particularly pronounced? 31 In the following we will attempt to answer this question.In doing so, we will explore the difference between MASH and FSSH in terms of the language of decoherence, revisiting the reasons for the breakdown of FSSH in systems with weak diabatic couplings, and showing how MASH improves upon these issues.We will begin by giving an overview of the two methods, highlighting the key similarities and differences between the FSSH and MASH algorithms.We will then describe how to simulate nonadiabatic rates using these approaches, before a detailed discussion of how each of the methods performs for a range of different physically relevant parameter regimes.
Methods: Here we give a brief description of the two methods in the case of a two-level system.Both FSSH and MASH treat the nuclear motion classically, with the nuclear positions and momenta represented by the classical variables q(t) and p(t) respectively.Between hopping events, the nuclei evolve under a force which is given by the derivative of the adiabatic potential corresponding to the 'active surface' where n is the active-state variable, and we label the upper adiabat + and the lower adiabat −.Electronic wavefunction coefficients, c±(t), are then propagated according to the time-dependent Schrödinger equation under the Hamiltonian generated by the nuclear trajectory.In both theories, these coefficients are used to determine when to hop, but are not used to calculate adiabatic population observables, which are instead obtained directly from the fraction of trajectories on a given active surface. 4An intuitive picture of the electronic dynamics can be obtained using the coordinates of the Bloch sphere This highlights the equivalence of the electronic dynamics to the rotation of a classical spin around a magnetic field where V± are the potentials corresponding to adiabatic states ϕ±, and dµ = ϕ+ ∂ϕ − ∂qµ is the nonadiabatic coupling vector.
What differs between FSSH and MASH is how the hops between surfaces are determined.Within FSSH, the probability of hopping from one surface to the other in a time-step δt is given by where negative probabilities indicate no hop.In contrast to this, the active surface in MASH is obtained deterministically by the simple condition i.e. the active state is the one with the larger probability, |c±(t)| 2 .The fact that MASH is deterministic might seem surprising, particularly given that it is the stochastic nature of FSSH that allows it to describe wavepacket splitting.However, as in other mapping-based methods, 39 the stochastic nature of surface hopping is replaced in MASH by sampling over initial wavefunction coefficients, as we shall explain below.To complete the specification of the dynamics, we need to define what happens to the momentum at a hopping (or attempted hopping) event.While there has been some debate in the literature as to how this should be done in FSSH, 47,48 the derivation of MASH from the QCLE leads to a unique prescription for how to deal with momentum rescaling and so-called frustrated hops (where the trajectory has insufficient energy to hop).The result is equivalent to what was originally argued for by Tully 49 (along with many others 50,51 ): the momenta are rescaled along the direction of the nonadiabatic coupling and are reflected in all cases that they do not have sufficient energy to hop.This suffices to describe the dynamical evolution of MASH and FSSH, however there is one additional important difference: how the simulation is initialized.For ease of comparison between FSSH and MASH, we will focus here on the calculation of correlation functions that involve only adiabatic populations and nuclear configurations (although we note that the MASH derivation leads to a rigorous prescription for the calculation of correlation functions involving electronic coherences).For a system starting in a specific adiabatic state, both FSSH and MASH are initialized with the corresponding active state, n(0).In FSSH the wavefunction coefficients are initialized as the corresponding pure state, e.g. if the initial state is n = + then c+(0) = 1 and c−(0) = 0 and the initial S vector points to the north pole of the Bloch sphere.In contrast the wavefunction coefficients in MASH are sampled such that the initial S is distributed over the entire hemispherical surface of the Bloch sphere corresponding to the initial state, with a probability density proportional to |Sz|.It is this sampling that effectively replaces the stochastic nature of the hops in FSSH.
Rate Calculations: Full details of the calculation of rate constants with MASH and FSSH are discussed in the supporting information.Here we give an overview of the most important aspects of reaction rate theory, focusing on the advantages of MASH over FSSH in two key areas: efficiency and accuracy.
Typically, the accurate determination of rate constants from a direct simulation of the population dynamics is not possible, as the barrier crossing is a rare event and prohibitively long trajectories would be required to observe a statistically significant number of reactions.The standard approach used to overcome this problem is the fluxcorrelation formalism. 52This avoids the rare-event problem by reformulating the rate in terms of a correction to transition state theory: the transmission coefficient.Importantly, the calculation of the transmission coefficient only involves running a short simulation up to the 'plateau' time, t pl , which is much shorter than the timescale of the reaction, t pl ≪ τrxn, but long enough that the initial transient behavior has subsided and the population decay is exponential. 52nfortunately, the FSSH dynamics do not obey timetranslation symmetry, and hence the flux-correlation formalism does not rigorously give the same result as calculating the rate from direct population dynamics.A number of approaches to overcome this issue have been suggested, such as using initial wavefunction amplitudes in the fluxcorrelation function generated from approximate backwardspropagation schemes, 34,53 as well as the use of dynamical enhanced sampling in the form of forward-flux sampling. 54ere, to avoid making further approximations, we simply calculate the FSSH rate from direct population dynamics, which is achievable due to the low computational cost of the model employed.The calculation of reaction rates with MASH presents a significant advantage in this regard: the dynamics of MASH do rigorously obey time-translation symmetry.This means that all of the usual machinery of the flux-correlation formalism (such as the Bennett-Chandler method [55][56][57] ) can be used to improve the efficiency of rate calculations in a way that is rigorously equivalent to the rate that would be obtained (less efficiently) with a direct simulation of the population dynamics.
2][33][34][35] This error is known to occur in problems where the system passes through regions of strong nonadiabatic coupling (equivalent to weak diabatic coupling) multiple times, resulting in an active state that is inconsistent with the wavefunction coefficients.Importantly, the dynamics of MASH can never become inconsistent in the way they do in FSSH, as the active state is determined explicitly from the wavefunction coefficients.This means that one may expect the overcoherence error to be less significant in MASH than it is in FSSH.In order to assess this, we consider the MASH and FSSH dynamics in two different regimes.Firstly, we focus on how the error affects dynamics near the plateau time, t ∼ t pl .This is done by calculating rates from the slope of the product population ⟨Pp(t)⟩ after the initial transient behavior has subsided for a system initialized in the reactant well in a classical thermal distribution.Secondly, we consider the dynamics over the timescale of the reaction, t ∼ τrxn, by simulating the full population decay.
Model: In order to compare numerically the accuracy of MASH and FSSH for the calculation of nonadiabatic rates, we consider the prototypical model for electron trans-fer: the spin-boson model. 58For ease of interpretation, we will consider the Brownian-oscillator form of the spin-boson model, [59][60][61][62] which consists of a harmonic (mass-weighted) solvent polarisation coordinate, Q, and an Ohmic bath describing the effect of friction along Q, with spectral density J(ω) = γω.The diabatic potentials along the solvent polarisation coordinate are then the famous Marcus parabolas 8 where Λ is the Marcus reorganisation energy, ε is the reaction driving force, and Ω is the characteristic frequency of the parabola.The two diabatic states are coupled by a constant diabatic coupling ∆, and the resulting adiabatic potentials along the solvent polarisation coordinate are given by (7) In both FSSH and MASH simulations, the initial positions and momenta are sampled from the classical Boltzmann distributions (not Wigner functions), and the initial active state is chosen with the associated Boltzmann weighting.As the nuclei are classical, the coupling of the solvent polarisation coordinate, Q, to its environment can be implemented efficiently using a Langevin equation with the friction coefficient γ.][61][62] In the limit of weak diabatic coupling (∆ → 0), Marcus theory predicts that the rate to go from one well to the other is given by where β = 1/kBT is the inverse temperature.Importantly Marcus theory is exact for this model in the weak-coupling limit under the assumption that the nuclear motion can be treated classically, i.e. in the absence of nuclear quantum effects such as zero-point energy and tunneling.This makes Marcus theory a very useful benchmark for assessing the accuracy of FSSH and MASH, which also assume the nuclear motion can be treated classically.In order to assess their behavior for intermediate values of ∆, where Marcus theory is not applicable, numerically exact quantum-mechanical rates were calculated using the hierarchical equations of motion (HEOM). 63,64All HEOM calculations were performed using the HEOM-Lab code 65,66 following the method described in Refs.62 and 67.
For both MASH and FSSH, the long-time behavior of ⟨Pp(t)⟩ is independent of the precise definition of reactants and products.* However, the definition of reactants and products will affect its short-time behavior.The optimum choice for the calculation of rates is the one for which the dynamics of a system initialized in the reactants most quickly settles into exponential decay.Normally this is a purely practical matter, however, choosing a (near) optimal definition has an additional importance in the present study: it allows us to separate the short-and long-time errors.The definition we use is that everything on the lower adiabatic surface to the right of the diabatic crossing or on the upper adiabatic surface on the left of the diabatic crossing is the product, and vice versa for the reactant.Mathematically this corresponds to where h(x) is the Heaviside step function and Pr = 1 − Pp.This definition works well for all the cases considered in this work.We demonstrate numerically in the supporting information that this gives the same rate constants as a purely position-space definition in the normal regime, or a purely adiabatic definition in the inverted regime, while having a shorter transient.
The parameters for the model are taken to be βΛ = 12, βℏΩ = 1/4, γ = Ω, for a range of values of ε and ∆.These parameters were chosen to allow a clear comparison of the accuracy of MASH and FSSH, at a reasonable computational cost.In particular the reorganisation energy was chosen to be high enough that the population transfer is in the slow incoherent limit, but low enough that it is possible to run direct population dynamics.This allows us to directly calculate FSSH rates, without needing to employ backward propagation or forward-flux sampling.Additionally, it allows us to demonstrate numerically that in MASH direct population dynamics are equivalent to the results obtained using the flux-correlation formulation, which we show in the supporting information.The characteristic frequency was chosen to make the system as classical as possible without the HEOM calculations becoming too expensive.This was done as our focus here is on assessing the relative accuracy of the dynamics of MASH and FSSH, rather than the importance of nuclear quantum effects.Finally it is known that the effect of overcoherence error becomes less pronounced at high friction, 35 and hence to make the test of MASH as stringent as possible we consider a system in the underdamped γ < 2Ω regime.Systems with larger reorganisation energy and higher friction are considered in the supporting information.
Results and Discussion: Figure 1 compares the rates calculated at the plateau time for a symmetric reaction, ε = 0, as a function of the diabatic coupling, ∆.We see that, for intermediate-to-large values of diabatic coupling, log 10 (β∆) ≳ −0.75, MASH, FSSH and HEOM all closely agree, with the HEOM rate showing only a slight ∼ 10% enhancement due to shallow tunneling.For smaller values of ∆, the reaction approaches the golden-rule regime where Marcus theory is valid.Here we see that MASH continues to closely match the exact results predicted by HEOM, while FSSH begins to deviate significantly with an unphysical slope.2][33][34] However, it raises the question: why does MASH not show a similar error?
To understand this, in Fig. 2 we analyse ⟨Pp(t)⟩ for the smallest value of ∆ considered in Fig. 1.The top left panel of Fig. 2 shows the full ⟨Pp(t)⟩.Although MASH and FSSH agree during the initial transient, the slope after this time differs significantly, with FSSH predicting a much slower population transfer.The remaining panels decompose ⟨Pp(t)⟩ into contributions from trajectories that have hopped 0 or 2 times between t = 0 and the current time, t.
The top right panel shows the sum of the zero-and two-hop trajectories.We see that the difference in the slopes of the MASH and FSSH curves closely resemble those in the full ⟨Pp(t)⟩, implying that other terms are contributing only to the transient and not the rate.Hence, to understand the difference between the MASH and FSSH rates, one can focus on just these trajectories.Unsurprisingly, the no-hop contribution to ⟨Pp(t)⟩ (which involves just a single passage through the crossing region) agrees very closely between MASH and FSSH.The key difference occurs in the trajectories that hop twice.The contribution of these trajectories, along with a depiction of a corresponding typical reactive path, is shown in the bottom right panel.From this we see that trajectories that hop twice contribute significantly (and correctly) to the rate in MASH but only contribute a very small amount in FSSH.Hence, the rate predicted by FSSH can be expected to be up to a factor of 2 too small, as previously pointed out by Jain and Subotnik in Ref. 34.
Having established that it is the two-hop trajectories that differ between MASH and FSSH, it remains to be explained why these trajectories go wrong in FSSH but not in MASH. Figure 3 illustrates the behavior of a typical two-hop trajectory in FSSH that "should" react but doesn't.The trajectory starts in the reactant well at t = 0.At t ≈ 7βℏ the trajectory reaches the crossing, and hops up due to the strong nonadiabatic coupling and correspondingly large hopping probability.Having hopped up, the trajectory then continues on the upper state before turning around and coming back towards the avoided crossing.Note, at this point the trajectory is not significantly affected by inconsistency or overcoherence error, as the wavefunction coefficients are essentially still in a pure state corresponding to the active sur-  face (i.e.Sz ≈ 1).At t ≈ 10βℏ the trajectory passes through the avoided crossing for a second time and most trajectories hop down (returning to the reactants).However, we follow one of the few that remain on the upper surface (probability ∝ ∆ 2 ).Now the wavefunction (which is predominantly in the lower state, Sz ≈ −1) is inconsistent with the active surface.When the trajectory returns to the avoided crossing for a third time, we expect it to hop down to the product well.However, the wavefunction is evolving in the opposite direction to the expected hop (from down to up instead of up to down).Hence, the probability to hop down is almost zero, and the trajectory incorrectly stays on the upper surface, leading to no reaction.In contrast MASH trajectories cannot have this problem.When an equivalent MASH trajectory approaches the avoided crossing for the third time, its spin vector is guaranteed to correctly point up (due to the consistency between its spin vector and the active surface).On passing through the crossing region, its spin vector will then flip down to the lower hemisphere, resulting in a downward hop and a successful reaction.So far we have only considered the dynamics on the timescale of a single barrier crossing.However, in the limit of weak diabatic coupling, the system may come back to the diabatic crossing (region of large nonadiabatic coupling) many times before the reaction takes place.This can lead to a build-up of overcoherence error, causing the long-time rate behavior to deviate significantly from the short-time behavior.To investigate this effect, Fig. 4 shows the population of products, ⟨Pp(t)⟩, for the full population decay, for two different driving forces, βε = 0 and βε = 3 = βΛ/4, with all other parameters kept the same as in Fig. 2.
Considering first the upper panel of Fig. 4, where βε = 0, we see immediately that the long-time behavior of both MASH and FSSH agrees perfectly with Marcus theory.This is a surprising result, as based on the short-time behavior we would expect FSSH to be too slow.However, it can be ex- The problem for FSSH occurs on the third crossing, where the wavefunction is inconsistent with the active state.Sz is then predominantly moving up, meaning that the probability to hop down is almost zero.This example is taken from a calculation with βℏΩ = 1/4, γ = 0, βε = 0, βΛ = 12.
plained away as a fortuitous cancellation of errors due to the symmetry of the model when ε = 0.This assumption is confirmed by considering the behavior of an asymmetric reaction, βε = 3, as shown in the lower panel.The short-time behavior of the symmetric and asymmetric systems are similar as can be seen from the inset.† At long time, however, we see that for the asymmetric system there is no fortuitous cancellation of errors.Instead, the build up of overcoherence error in FSSH leads to a population decay that is noticeably too fast, with a half-life approximately 3.5 times shorter than Marcus theory.In contrast, MASH goes from being almost exact at short time to being about 1.4 times too fast at long time.‡ We see, therefore, that while both MASH and FSSH suffer from a build-up of overcoherence error at long time this error is significantly more pronounced in FSSH.
The build-up of overcoherence error at long time is of course well established.While it is nice that this error is much smaller in MASH than FSSH, in real simulations on such incredibly long timescales one should apply decoherence corrections in both theories.In this regard, the short-time accuracy of MASH also presents a significant advantage.To see why, we note that the application of decoherence corrections is always a balancing act: you need to apply them often enough to fix the overcoherence error, but apply them too often and you will force the system to remain forever on the same adiabat (the quantum Zeno effect).The advantage  of MASH is that it requires decoherence corrections less often to obtain accurate results.This means that they can be applied only in regions where it is safe to do so, such as the reactant wells, and not in the vicinity of the coupling region.This makes it more robust and means simpler decoherence schemes can be successfully used.
In Fig. 4 we demonstrate this by considering the behavior of MASH and FSSH when one applies a simple decoherence correction.(Note decoherence corrections in MASH correspond to resampling S from the hemisphere corresponding to the current active state according to the |Sz| weighting. 36) For simplicity we use an energy cutoff such that decoherence corrections are applied only when the gap between the states is large, V+ − V− > 4kBT , i.e.where the system is far from the avoided crossing.§ In the upper panel we see that, for the symmetric system ε = 0, the population decay predicted by MASH is unaffected by application of the decoherence correction, leaving it in perfect agreement with Marcus theory.In contrast, however, the FSSH results are made significantly worse by application of the decoherence correction, for reasons explained below.For the asymmetric system, βε = 3, we see that application of the decoherence correction, improves the original MASH result, removing the ∼ 40% error, and bringing it into almost perfect agreement with Marcus theory.Again, however, the simple decoherence correction does not fix FSSH, in this case taking the rate from being too fast to too slow.These results can be understood by noting that by applying the decoherence correction far away from the crossing region we simply make the long-time dynamics consistent with the short-time dynamics.For MASH the short-time dynamics has the correct rate, but for FSSH the short-time dynamics is wrong, as can be seen from the inset, and hence we recover the spuriously low rate seen in Fig. 1.That such simple decoherence corrections do not fix FSSH is not a new observation, and for this reason many far more sophisticated decoherence approaches have been developed. 28,29However, these methods often come with additional disadvantages, such as increased cost, and as they are ad-hoc they are not always guaranteed to improve the results.The point we would like to stress here is that the increased accuracy of MASH at short times means that decoherence corrections can be applied much more infrequently.For many ultrafast problems this means that they may not be needed at all.But when they are needed they can be both safer and simpler.
Finally, having understood the difference between MASH and FSSH, we consider the famous Marcus turnover curve.Figure 5 shows the behavior of the rate in the weak-coupling limit (log 10 (β∆) = −7/5) as a function of the bias to products, ε.As in Fig. 1 the FSSH and MASH rates are calculated from the slope of ⟨Pp(t)⟩ near the plateau time, between t = 10βℏ and t = 20βℏ.As expected from the results above, FSSH deviates significantly from Marcus theory and the exact results, showing an unphysical asymmetry about ε = Λ.In contrast MASH reproduces both the exact results and Marcus theory very well.MASH is also not perfectly symmetric due to a slightly larger error deep in the inverted regime (ε = 2Λ) than in the symmetric case (ε = 0).However, in both cases the errors are less than 10%.The largest error in MASH is observed close to the activationless limit ε/Λ = 1.Here the MASH rate is about 15% higher than the Marcus-theory result.In contrast FSSH is about 60% too large.As the avoided crossing is located at the minimum of the reactant well in the activationless case, this leads to a faster build up of overcoherence error and the increase in the rate seen at long time in the lower panel of Fig. 4 starts to affect the dynamics even at the short times considered here.This is confirmed by application of the same decoherence correction as was used in Fig. 4, which stops the build up of overcoherence error, resulting in MASH rates that are within 7% of the exact rate for the full range of ε considered.As in Fig. 4 the increased error of FSSH at short time means that this simple decoherence correction is not sufficient to fix the inconsistency error of FSSH, and hence the rates still deviate significantly from the Marcus-theory result, as can be seen in Fig. S6 of the supporting information.
2][33][34][35] Here we have revisited this problem to assess the accuracy of a newly proposed alternative to FSSH, the mapping approach to surface hopping (MASH).In comparing MASH and FSSH, we have considered two different timescales: the timescale of a single barrier-crossing event, t pl , and the timescale of the reaction, τrxn.
On the timescale of barrier crossing, MASH provides a significant improvement upon FSSH, accurately recovering the results of Marcus theory without the use of decoherence corrections.This might seem surprising at first, as it is not immediately obvious how MASH, which is also an independent trajectory method, is able to capture decoherence.However, we have shown that the improvement can be explained in terms of the dramatic inconsistency between the active state and wavefunction coefficients, which can exist in FSSH but is absent from MASH.
On very long timescales, MASH again provides a significant improvement over FSSH.While overcoherence error does still build up in MASH, we have found it to be much less significant than in FSSH.This can again be explained in terms of the inconsistency in FSSH, which means that the build-up of error can be sudden and large, whereas in MASH the build-up of error is more gradual and ultimately smaller.Perhaps most importantly, the increased accuracy of MASH over FSSH at short time means that, when they are used, decoherence corrections need only be applied well away from the coupling region, making them safer and simpler to use.
MASH also comes with additional practical advantages over FSSH in the calculation of reaction rates.In particular, as the dynamics of MASH are deterministic and obey time-translation symmetry, there is no need for approximate backwards-time propagation or advanced methods such as forward-flux sampling.One can instead rigorously apply the flux-correlation formalism and related techniques, such as the Bennett-Chandler method, in order to efficiently calculate reaction rates.Given these significant improvements, and the fact that MASH is simple to use, requires only relatively minor modifications to existing FSSH code, and can be run at equivalent computational cost, MASH has the potential to replace FSSH as the go-to method for the simulation of nonadiabatic processes.
The only thing limiting MASH as a replacement to FSSH is that the current theory is restricted to two-state problems.Recently a modification to MASH has been proposed, designed for application to multistate problems. 70However this theory is a different method to the MASH described here.It does not reduce to the current theory in the case of a two-level system, and although it is accurate for many problems, it was shown to be significantly less accurate for the timescales of population decay in a spin-boson model in the inverted regime.Work to develop a multistate generalisation of the present MASH method is in progress, and if this can be achieved, while retaining the advantages of the two-state theory, it would present a significant challenge to the hegemony of FSSH.
Finally, we note that we have focused here exclusively on the limit of classical nuclei.0][81][82][83][84][85] While there has been significant development in methods specialised for accurately predicting thermal reaction rates, [86][87][88][89][90][91] at present there is no fully dynamical method that can offer comparable accuracy. 92This is in part due to the difficulty that such dynamical methods face in accurately describing rates even in the limit of classical nuclei.In this regard, the results of the present study indicate that MASH provides a new and exciting route to the development of a fully dynamical nonadiabatic theory capable of accurately describing nuclear tunneling and zeropoint energy.
Supporting Information: Details on rate calculations and additional results to support the conclusions of the paper.Secondly, in the inverted regime we compare to using a purely adiabatic state definition of the reactants and products Figure S1 compares ⟨P p (t)⟩ calculated using the general definition (used in the main paper) against the position and adiabatic definitions.The lower panel shows the results for a system in the normal regime and the upper panel shows the results for a system in the inverted regime.We see in both cases that at long time the definition used in the main paper is equivalent to the more commonly used definitions given in Eqs.S14 and S15.
of Eq.S27 corresponds to whether a hop is attempted and the second to whether it has enough energy to actually occur.The first term is only non-zero if the current rate of change predicts a change in the sign of S z between t and t+ϵ (i.e. an attempted hop).The second term then accounts for the possibility that the hop is frustrated, in which case the overall expression must be zero (due to the reflection of S z off the spin-sphere equator for a frustrated hop).To see that this is the correct expression, note that if the system is on the upper surface, h(S z (t)) = 1, then the hop is necessarily not frustrated (this is true whether t is just before or just after a hop, i.e. whether ϵ is positive or negative).However, if the system is on the lower surface, h(−S z (t)) = 1, then it will hop (or just have hopped) if and only if the kinetic energy along the derivative coupling vector is greater than the adiabatic energy gap, ∆E hop (t) > 0.
As we have already alluded to, there is a choice in whether the limit is taken from above (ϵ → 0 + ) or from below (ϵ → 0 − ).First we will take the limit as ϵ → 0 − .This will result in expressions for which t is infinitesimally after the hopping time, rather than before.Noting that only the first term on the right-hand side of Eq.S27 involves ϵ, the key to evaluating the limit is the following identity lim where we will replace x and b with S z (t) and Ṡz (t).This result is straightforward to verify by considering the two cases b > 0 and b < 0 separately.Combining this with Eq.S27, and utilising trivial identities such as δ + (S z (t))h(−S z (t)) = 0, one then obtains the following expression for the time derivative This expression has the advantage that it is symmetric, although numerical implementation will involve trajectories that hop during the first time step, which one might want to avoid.Alternatively, therefore, one can make use of Eq.S31, which has the hop infinitesimally after t.This might seem like it is the wrong choice, but after inserting into Eq.S20 and making use of the time-inversion symmetry the hop will be infinitesimally before t = 0, resulting in the following expression  Numerical implementation of flux-correlation formalism.

Position based definition (normal regime).
With some straightforward algebraic manipulations, Eq.S23 can be decomposed in the usual Bennett-Chandler form S6 as follows where ⟨. . .⟩ Q(0)=Q ‡ indicates an average over the distribution Note that in this form, it is assumed that S z is sampled with the correct Boltzmann weights.
Alternatively, we can sample S uniformly from the Bloch sphere and write where ⟨. . .⟩ Q(0)=Q ‡ ,MF indicates an average over the distribution Adiabatic state defintion (inverted regime).
We can evaluate Eq.S33 using a modified Bennett-Chandler procedure.We start by noting that the presence of the derivative coupling in Ṡz = − j 2d j (q)p j m j S x (S44) means that trajectories starting in regions of high derivative coupling will dominate.Ideally, one should therefore incorporate this into the sampling.Practically for our purposes, the sampling of the nuclear positions and momenta can be done from an as yet undefined distribution ρ samp (p, q) = e −βHsamp(p,q) tr e −βHsamp(p,q) .(S45) The first term of Eq.S33 can be written in terms of averages over this distribution as The integrals over the spin-sphere in the second term can be performed analytically to give such that overall we have

S14
For the term constrained initially to the lower adiabatic surface, we can follow the same steps to obtain the final equation for the derivative of the population as Here we give the expression for the symmetric definition of the derivative, Ṗ (a) p , however a similar expression can also be derived for Eq.S34.Note that sampling from the distribution with S z (0) = 0 + or S z (0) = 0 − is implemented in practice by just setting the initial S z to a small floating point number, ±10 −10 .

General definition (used in main paper).
Combining the ideas of the last two sections, we can arrive at a Bennett-Chandler scheme to efficiently compute the rate with the general dividing surface.For the first term in Eq.S38, this can be achieved by writing Alternatively sampling S uniformly from the surface of a sphere, we have that For the second term and third term, we can follow the same approach as with the pure adiabatic S15 dividing surface.This gives Here we have chosen to use the symmetric definition of d dt h(S z (t)) (Eq.S32).This means that hops can occur within the first time step, however we found that this did not present a significant numerical difficulty.We found that it was more important to sample the initial momenta such that for each p, we also sampled a trajectory with −p at the same S and q.We note that this could have also been achieved by using the asymmetric version, where hops are always infinitesimally before t = 0, by appropriately rescaling the momenta.Finally in the present calculations, we used a very simple choice for H samp corresponding to sampling from the reactant diabatic potential origin, shifted to the diabatic crossing seam H samp (p, q) = f j=1 p 2 j 2m j + U 0 (q − q shift ) (S54) with q shift = (Q ‡ , 0, 0, 0, . . .).This is not the optimal choice, as it does not take into account the narrowing of the derivative coupling as ∆ → 0, however it was sufficient for the present purpose.
Additional Results.
In the following, we include results of additional systems, that support the conclusions of the main text. .FSSH and MASH rates were calculated from the slope of ⟨P p (t)⟩, at the plateau time, between t = 10βh and t = 20βh.This shows essentially the same behaviour as seen for the symmetric system in Fig. 1 of the main text.

S16
In Fig. S4, we show that the qualitative behaviour is unchanged from Fig. 1 of the main text when one adds a small bias to products, βε = 3.
In Fig. S5, we consider a system which is equivalent to that considered in Fig. 1 of the main text, except that it is in the overdamped γ = 4Ω rather than underdamped regime.This illustrates that the overcoherence problem is less significant at high friction.
In Fig. S6, we show the results of including the simple decoherence correction on FSSH.These were left out of the main paper to avoid clutter.This illustrates that unlike for MASH, the simple decoherence correction is unable to fix FSSH.
In Fig. S7, we consider a system with a large reorganisation energy, βΛ = 60, and high friction, γ = 32Ω.This illustrates the utility of the flux-correlation formalism, as direct calculation of the rates would be prohibitively expensive at the smallest values of ∆ considered.

Figure 3 :
Figure 3: An example of a typical incorrect "two hop" FSSH trajectory that fails to react, along with a comparable but correct MASH trajectory.The problem for FSSH occurs on the third crossing, where the wavefunction is inconsistent with the active state.Sz is then predominantly moving up, meaning that the probability to hop down is almost zero.This example is taken from a calculation with βℏΩ = 1/4, γ = 0, βε = 0, βΛ = 12.

Figure 4 :
Figure4: Full population decay for a symmetric, βε = 0, and asymmetric, βε = 3, spin-boson model, with βℏΩ = 1/4, γ = Ω, βΛ = 12 in the limit of weak diabatic coupling, log 10 (β∆) = −7/5.Inset shows ⟨Pp(t)⟩ at short to intermediate time.Decoherence corrections are applied only when the energy gap is large (V+ −V− > 4kBT ), making long-time behavior consistent with short to intermediate time.This illustrates that not only is MASH more accurate than FSSH without the application of decoherence corrections, but also that, unlike in FSSH, simple decoherence corrections are sufficient to bring MASH into line with the correct result.

Figure S2 :
Figure S2: Demonstration of the MASH flux-correlation formalism.Upper panel: comparison of direct population dynamics (Eq.S11) and the cumulative integral of the flux-correlation formalism (Eq.S38) for the symmetric, βε = 0, spin-boson model with log 10 (β∆) = −1, βhΩ = 1/4, γ = Ω, and βΛ = 12.Lower panel: generalised flux-correlation function calculated using the fluxcorrelation formalism (Eq.S38) the long-time limit of which defines the rate constant.All results use the same spin-boson model, and the general definition of reactants and products.

Figure S4 :
FigureS4: Log-log plot of the rate against the diabatic coupling, for an asymmetric, βε = 3, spinboson model, with βhΩ = 1/4, γ = Ω and βΛ = 12.FSSH and MASH rates were calculated from the slope of ⟨P p (t)⟩, at the plateau time, between t = 10βh and t = 20βh.This shows essentially the same behaviour as seen for the symmetric system in Fig.1of the main text.