A multi-state mapping approach to surface hopping

We describe a multiple electronic state adaptation of the mapping approach to surface hopping introduced recently by Mannouch and Richardson (J. Chem. Phys. 158, 104111 (2023)). This adaptation treats populations and coherences on an equal footing and is guaranteed to give populations in any electronic basis that tend to the correct quantum-classical equilibrium values in the long-time limit (assuming ergodicity). We demonstrate its accuracy by comparison with exact benchmark results for three- and seven-state models of the Fenna-Matthews-Olson complex, obtaining electronic populations and coherences that are significantly more accurate than those of fewest switches surface hopping and at least as good as those of any other semiclassical method we are aware of. Since these results were obtained by adapting the scheme of Mannouch and Richardson, we go on to compare our results with theirs for a variety of problems with two electronic states. We find that their method is sometimes more accurate, and especially so in the Marcus inverted regime. However, in other situations the accuracies are comparable, and since our scheme can be used with multiple electronic states it can be applied to a wider variety of systems.


I. INTRODUCTION
Fewest switches surface hopping (FSSH) is the dominant method for simulating electronically nonadiabatic dynamics in chemistry.3][4][5] It remains popular for three (good) reasons: (i) it is easy to use; (ii) it is robust when used with ab initio potentials; and (iii) it has been found to work reasonably well in many applications.However, while there have been attempts to derive it from first principles, 6 FSSH is still widely considered to be an ad hoc algorithm.It also suffers from a well-known 'overcoherence' problem that has led to the development of a number of different 'decoherence' corrections, 7 none of which has become universally accepted.
A second school of thought advocates treating the electronic and nuclear degrees of freedom on an equal footing by mappling the multi-state electronic system onto a classical phase space.The simplest approach of this type is (multiple-trajectory) Ehrenfest dynamics, which is based on a mean-field coupling of the nuclear and electronic motions.However, this is known to produce various unphysical results including an overheating of the electronic subsystem. 8Semiclassical mapping methods represent the multiple electronic states as harmonic oscillators 9,10 or generalized spins, 11 and use either initial phase space sampling or symmetrical quasi-classical binning 12,13 to calculate the dynamics.These mapping methods typically give better accuracy than either FSSH or Ehrenfest dynamics in applications to system-bath models of condensed-phase systems.However, they do so at the cost of allowing the electronic-state populations to a) Electronic mail: johan.runeson@chem.ox.ac.uk become negative, which can result in the nuclei evolving (and potentially even diverging) on inverted potentials.While this problem may not be visible in simple models, it does help to explain why mapping methods are only rarely combined with ab initio potentials and applied to realistic chemical systems.
Given the complementary strengths and weaknesses of surface hopping and phase-space mapping, it is natural to ask whether they could be combined to give a better method.Very recently, Mannouch and Richardson have explored this idea and used it to develop a new 'mapping approach to surface hopping' (MASH). 14As in surface hopping, this method uses trajectories that hop between the physical adiabats, thereby eliminating concerns about dynamics on inverted mapped potentials.But rather than employing stochastic hops, the active surface is obtained deterministically from the electronic wavefunction, as it is in the phase-space mapping approach.This has several distinct advantages, not least of which is that it avoids the need for heuristic 'decoherence' corrections.The results can instead be improved (if necessary) through a careful resampling of the electronic wavefunction (the quantum-jump procedure 15 ).Even without this resampling, MASH has been found to give better results than either pure surface hopping or pure mapping for a wide range of two-state problems, with no more computational effort. 14It is thus arguably the most promising method that has yet been proposed for nonadiabatic dynamics.
However, the formulation of MASH that Mannouch and Richardson presented was restricted to two coupled electronic states. 14It is not obvious how to extend it to more states because they based their formulation on correlation functions involving specific two-state prescriptions for the electronic populations and coherences that were specialized to the adiabatic basis (see Sec II).A general nonadiabatic dynamics method should be applicable to an arbitrary number of electronic states, it should treat populations and coherences on an equal footing, it should provide results that transform correctly under the unitary rotations that take one electronic basis to another, and it should be able to provide these results directly in any chosen basis.In Sec.III, we describe such a multi-state adaptation of MASH, and show that it is guaranteed to give the correct quantum-classical equilibrium populations of the electronic states in any basis in the long-time limit (assuming ergodicity).In Sec.IV, we demonstrate that this adaptation works just as well for a standard multi-state exciton energy transfer problem as Mannouch and Richardson have shown MASH to work for two-state problems. 14In Sec.V we compare our adaptation with Mannouch and Richardson's original version of MASH for a variety of problems with just two electronic states, and in Sec.VI we conclude this paper.

II. NONADIABATIC DYNAMICS
Consider a general nonadiabatic system defined by the Hamiltonian where q = {q j } and p = {p j } are the coordinates and momenta of f nuclear degrees of freedom with masses {m j }, and V (q) is the potential energy operator.In model systems, it is usually most convenient to express the potential in a given diabatic (q-independent) basis {|n⟩}, i.e., V (q) = nm V nm (q)|n⟩⟨m|. ( Although it would in principle be possible to evolve surface hopping dynamics directly in the diabatic basis, it is in practice almost always run in the adiabatic basis (the local eigenbasis of V (q)), where The two bases are related through their transformation matrix elements U na (q) = ⟨n|a(q)⟩.In this article, we will consider model systems defined in a diabatic basis, but our method can also be used directly with potentials in the adiabatic basis, such as those provided by ab initio electronic structure calculations.For most photochemical applications, it is a reasonable approximation to treat the nuclei as classical particles.However, to consistently evolve a coupled system of classical and quantum degrees of freedom is a longstanding problem in semiclassical dynamics.To see why, consider a particle at configuration q with momentum p and with electronic state |ψ⟩ = n c n |n⟩.The natural starting point is to evolve the electronic state according to Schrödinger's equation of motion.This can be done equivalently either in the diabatic basis or in the adiabatic basis where d j ab (q) = ⟨a(q)|∇ j |b(q)⟩ is an element of the nonadiabatic coupling vector.(Throughout this paper we use units where ℏ = 1.)The main difficulty arises when constructing the nuclear dynamics, and in particular when considering the 'back-action' of the electrons on the nuclei.There is a 'force operator' Fj (q) = −∇ j V (q), but this is not yet useful to run classical dynamics.What one would like is equations of motion of the form where F j (q) is a (yet to be defined) classical force.Existing schemes for nonadiabatic dynamics differ mainly in the way they construct this force.In the following subsections we briefly summarize and comment on some of the most important strategies.

A. Fewest switches surface hopping
In surface hopping, the instantaneous force on the nuclei is taken to be that of a single adiabatic state, called the active surface.If the active surface is a, then this force is F j (q) = −⟨a(q)|∇ j V (q)|a(q)⟩.(7)   If the trajectories enter a region with non-zero nonadiabatic coupling, they can switch active surface (or 'hop').
At each discrete time step, the probability to hop from a to b is taken to be Tully chose this probability such that the fraction of trajectories evolving on surface a would approximate the average population ⟨|c a | 2 ⟩ with a minimal number of switches. 1If a hop occurs, the momentum is rescaled along the nonadiabatic coupling vector such that the total energy is conserved.If there is not sufficient kinetic energy to overcome the difference in potential energy between the pre-and post-hop adiabatic surfaces, the standard (although not universally accepted 16 ) practice is to reject the hop and reverse the momentum in the direction of the nonadiabatic coupling vector.Despite the considerations behind the choice of hopping probability, the fraction of trajectories evolving on surface a does not strictly agree with ⟨|c a | 2 ⟩.This inconsistency is at the root of many of the issues present in surface hopping. 17,18Traditionally, the problem has been identified as an 'overcoherence' of the electronic coefficients that can be overcome with (more or less heuristic) 'decoherence corrections'.Many such corrections have been proposed, but despite much effort, there is as yet no consensus as to whether any of them has solved the underlying problem.Furthermore, surface hopping does not generally guarantee relaxation to the correct longtime equilibrium in condensed-phase systems, although it does often provide a better approximation than Ehrenfest dynamics. 19,20

B. Semiclassical mapping approaches
A rather different strategy is to construct the force to be a coherent average over contributions from multiple electronic states.A simple approach of this type is Ehrenfest dynamics, which uses the expectation value of the force operator, Ehrenfest dynamics has several severe drawbacks, of which the most important are that it violates detailed balance (it relaxes the system to an overheated equilibrium) and fails to capture wavepacket branching in scattering models.Nevertheless, it has the advantage of being invariant to a unitary transformation of the electronic basis, and the deterministic nature of the force allows an ergodic analysis of its long-time limit.Explicitly, the real and imaginary parts of the electronic coefficients can be regarded as phase-space variables on the same footing as the nuclear degrees of freedom.In terms of these variables, it is clear that Eq. ( 4) is equivalent to the dynamics of a set of N harmonic oscillators.The last of these observations has inspired a more formal mapping of the electronic states to quantum harmonic oscillators, which is now known as the Meyer-Miller-Stock-Thoss (MMST) mapping. 9,10Taking the classical limit of this oscillator model leads to a classical phase-space theory with its own force as well as expressions for population and coherence estimators.The new force is also of coherent-average type, but differs from the Ehrenfest force in that the instantaneous populations can be negative (or larger than one).Despite this seemingly unphysical behaviour, the weighted average over many trajectories has for many model problems been found to be more accurate with the MMST mapping than in (multi-trajectory) Ehrenfest dynamics. 13,21armonic oscillators are not the only way to map the electronic coefficients onto classical variables.For twolevel systems, another choice would be to use the wellknown isomorphism to a spin-1/2 system.This approach was recently used to develop a 'spin mapping' analogous to the MMST mapping, 22 which leads to a subtly different definition of the force and the estimators.For N -level systems, the spin mapping has a natural generalization in terms of the so-called Stratonovich-Weyl transformation. 11This has been shown to at least partly solve the overheating problem of Ehrenfest dynamics, in the sense that the long-time equilibrium reduces to phase-space averages that agree with quantum mechanics up to first order in β = 1/k B T . 23owever, just like the MMST mapping, spin mapping can predict (unphysical) negative populations at low temperatures.For an individual trajectory, an instantaneous negative population can cause it to evolve and diverge on an inverted potential, which is completely unphysical.A recent attempt to overcome this problem was to construct a mapping to an anisotropic spin. 24This does remove the issue of inverted potentials and can (at least for two-level systems) be constructed so as to give the correct long-time equilibrium populations at any temperature, as long as the classical nuclear assumption is valid.However, the timescale of the relaxation to equilibrium was found to be worse in several cases than that of the original spin mapping.

C. Mapping approach to surface hopping
It should be clear from what we have said so far that the surface hopping and mapping approaches each have their own advantages and disadvantages.Given this, the natural question is whether it is possible to develop a method that combines the strengths and eliminates the weaknesses of the two strategies.
Mannouch and Richardson have recently shown that, at least for two-level systems, this may indeed be possible. 14Their 'mapping approach to surface hopping' (MASH) uses the nuclear force of a single active adiabatic surface as in surface hopping, but instead of treating the active surface V a (q) separately from the electronic wavefunction, they set it to be that of the adiabatic state with the largest instantaneous population |c a | 2 .In this way, there is no need to introduce a stochastic hopping probability, because the active state is always uniquely determined by the electronic wavefunction.
Figure 1 illustrates the situation on the Bloch sphere of a two-level system.Adiabatic states 1 and 2 correspond to the opposite poles along the axis parallel to V = (V x , V y , V z ), where V α = Tr[ V σα ] and the {σ α } are the Pauli spin operators.In MASH, the instantaneous active surface is set to 1 when the Bloch vector σ = ⟨ψ| σ|ψ⟩ is on the hemisphere closest to adiabatic state 1, and 2 when it is on the hemisphere closest to adiabatic state 2.
In addition to providing a deterministic alternative to stochastic surface hopping, the MASH approach has several other appealing features, a more detailed discussion of which can be found in Ref. 14: (i) it satisfies detailed balance 24 in the sense that it gives the correct long-time populations of the adiabatic electronic states (assuming ergodicity); (ii) the momentum rescaling and momentum reversal arise naturally in the deterministic surface hopping algorithm without any ambiguity; Any selection of two opposite fixed points on the sphere (e.g., the black dots) defines a diabatic basis, whereas the adiabatic basis corresponds to the opposite points on an axis that follows (Vx, Vy, Vz) (blue and red dots).In MASH, the nuclei evolve on the adiabatic state that is closest to the instantaneous direction of the Bloch vector, σ = ⟨ψ| σ|ψ⟩.
(iii) in place of heuristic decoherence corrections, there is a way to systematically improve the results by carefully resampling the electronic wavefunction along the trajectory (the quantum jump prodecure).
Mannouch and Richardson have demonstrated by comparison with quantum benchmark calculations that MASH is more accurate for typical system-bath models than both FSSH and state-of-the-art mapping approaches.Their method can treat wavepacket branching just at least as well as FSSH, which is currently only possible in mapping via a more expensive cancellation of positive and negative phase-space contributions. 25They also found MASH to be the most accurate classical-trajectory method considered so far in describing ultrafast internal conversion in pyrazine. 14owever, their formulation of the method was restricted to two-level systems, 14 and it is not obvious how to generalize it to more electronic states because Mannouch and Richardson specifically constructed their observables for the case of two.Explicitly, for electronic observables Â and B, they computed the time-correlation function 26 where ⟨•⟩ is an expectation value taken with respect to an appropriate density function over (p 0 , q 0 , c 0 ).Here the estimators A(c) and B(c) are constructed from the following complete set of mappings for two-level populations and coherences defined in the adiabatic representation, Although this scheme has been shown to work well for a multitude of problems, it is worth asking whether there exists a simpler procedure, especially if one aims to generalize it to more than two electronic states.Firstly, in our view, there should be no reason to give coherences a special treatment.Since these are just population differences in a rotated basis, one should be able to treat them on an equal footing to the populations.Secondly, it would be preferable to calculate the observables in any diabatic basis directly from the wavefunction coefficients in that basis, without having to first convert them to a linear combination of adiabatic observables.The adaptation of MASH presented in the next section achieves both of these goals while at the same time generalising the method to an arbitrary number of coupled electronic states.

III. MULTI-STATE MAPPING
Having set the context of what we are about to do, we are now ready to propose a multi-state generalisation of MASH.We will separately address the issues of how to run the dynamics, how to measure observables, and how to set up the initial conditions for a typical photochemically initiated nonadiabatic problem.

A. Dynamics
The dynamics extends quite naturally from the twolevel case by using the same expression for the nuclear force as in Eq. ( 7) and picking the active surface to be that of the adiabatic state with the largest instantaneous population.To express this more precisely, we define the N populations and phases in a given basis as If we also define the classical state projectors then precisely one of these N projectors will be non-zero in any basis at each point in the classical c space.As in two-level MASH, we take the dynamically active potential energy surface -the surface that is used to define the nuclear force in Eq. ( 7) -to be that of the unique FIG. 2. Schematic depiction of a MASH trajectory for two (left) and three (right) states.In any given basis, the set of N populations is a point on a simplex (generalized triangle) with the N unit vectors as vertices.In MASH, the active surface is set to be the adiabatic state with the highest instantaneous population.Thereby, each adiabatic state corresponds to the region on the simplex closest to a given vertex.Hops occur whenever the system crosses a border between two of these regions.
adiabatic basis state a with a non-zero classical state projector Θ a .This choice of the adiabatic potential for the nuclear dynamics is consistent with experience from decades of successful surface-hopping calculations.
Assuming the initial electronic wavefunction is normalized, the populations will always sum up to one in any basis, and they will always be non-negative, P n ≥ 0. (The electronic evolution in Eq. ( 4) is unitary and preserves these two conditions for all time.)In the two-level case, one can visualize (P 1 , P 2 ) as a point evolving on the line segment between (1, 0) and (0, 1).Similarly, in the three-level case (P 1 , P 2 , P 3 ) is a point on the triangle with vertices (1, 0, 0), (0, 1, 0), and (0, 0, 1). Figure 2 illustrates the regions of different active surfaces for these two cases.In general, for N levels, (P 1 , P 2 , . . ., P N ) is a point on the simplex which is a geometrical object that can be thought of as a high-dimensional triangle with vertices on the unit vectors. 27hen a trajectory passes from the region associated with adiabatic state a to that associated with state b, there is a 'hop', and we rescale the momentum along a direction determined by the nonadiabatic coupling vector as explained in detail in Appendix E. If the kinetic energy is insufficient to hop, we reverse the momentum along the same direction.So far, everything is just a direct extension of the two-level case. 14In the three-level case, however, we need to address a new issue, which is that hops can occur between states that are uncoupled from each other.For example, if states 1 and 2 are coupled but both are uncoupled from state 3, then trajectories that start in the blue region but are close enough to the green can hop to state 3 before they reach state 2. It is not obvious how to deal with this situation.The simplest option is to accept it.Another option would be to reject the hop if the coupling between the two relevant states (in this case 1 and 3) is less than a given threshold.As in the case of insufficient kinetic energy, one would in such a situation also reverse the momentum.Here we will adopt the simplest option and accept the hop regardless of how strongly coupled the states are.This avoids the ambiguity of defining a coupling threshold and it finishes our discussion of the dynamics.What remains to be decided is how to measure observables.Since our aim is a method that gives the correct quantum-classical equilibrium populations in the long-time limit, we shall begin by discussing what happens when the dynamics has reached equilibrium.

B. Equilibrium
Because MASH is deterministic (in contrast to standard surface hopping), one can make more powerful statements about its long-time limit.Provided the system is sufficiently ergodic (as is typically true in the condensed phase), any initial distribution will in the longtime limit reduce to the equilibrium distribution ρ(p, q, c) ∝ e −βE(p,q,c) , (15)   where the energy function is with T (p) = j p 2 j /2m j and Θ a ≡ Θ a (c).This means that one can predict the long-time values of all observables in terms of phase-space averages with respect to this distribution.Such an analysis has recently been done in the two-state case for a variety of mapping methods, including MASH. 28In the following, we consider equilibrium expectation values in MASH for a general N -level system.
For a coupled electronic-nuclear system, the quantummechanical partition function with classical nuclei is This expression mixes a classical phase-space integral with a quantum trace.MASH (like other mapping approaches) effectively replaces the trace with an integral over all normalized electronic wavefunctions c, where N = 2π N /N !normalizes the electronic integral such that |c|=1 dc/N = N (see Appendix A).Because each Θ a is an idempotent projector onto a different region of c space, and it has a unit integral |c|=1 dc/N Θ a = 1, it follows that which shows that MASH is consistent with the mixed quantum-classical partition function (Z MASH = Z).By the same argument, the thermal equilibrium population of each adiabatic quantum state, will be equal to the corresponding MASH expression provided Θ a is a state projection in the adiabatic representation.Hence, if one were to choose to use the state projectors as population estimators, this would be guaranteed to give the correct equilibrium populations in the adiabatic basis.However, it would not work more generally, because ⟨Θ n ⟩ does not equal ⟨P n ⟩ in other bases.To see this, we can write the latter explicitly as whereas the former is This is in general not equal to In other words, the state projection defined by Eq. ( 13) does not transform correctly as a population estimator under unitary basis rotations.

C. Populations
A population estimator Φ n can in fact be constructed such that ⟨Φ n ⟩ = ⟨P n ⟩ in any basis of N orthonormal electronic states, whether they are adiabatic or diabatic.
Here we shall give an overview of the argument that establishes this and leave the technical details to Appendices B and C.
Since the Ehrenfest populations P n = |c n | 2 do transform correctly (i.e., quantum mechanically) under unitary basis rotations, as do scalars, our ansatz is an equivariant population estimator of the form in which α N and β N are scalar parameters that are to yet be determined.Our goal is to choose these parameters such that holds for all adiabatic states a and diabatic states n.This will suffice to ensure -by virtue of the argument in Eqs. ( 22) and ( 23) -that the expectation value of Φ n will be correct at equilibrium.Appendix B shows that Eq. ( 26) will be satisfied if the simpler condition holds for all adiabatic states a and b.The constraints that this imposes on α N and β N are then investigated in Appendix C, which obtains the unique solution where H N = N n=1 1/n.Hence we arrive at the estimator which has the pleasingly democratic interpretation of measuring the population of state n relative to the situation in which all the populations are equal to 1/N (as they are for the states at the centre of the simplex).

D. Coherences
In contrast to Ref. 14, as well as most other literature on surface hopping, we shall begin by rewriting coherences as population differences between rotated states so that they can be treated in the same way as populations.Note that we are in an ideal position to do this here because our population estimator in Eq. ( 29) transforms correctly under unitary basis rotations.
The coherences between states |n⟩ and |m⟩ are linear combinations of the Hermitian operators σx nm = |n⟩⟨m|+ |m⟩⟨n| and σy and we can rewrite σx nm and σy nm as differences between the population operators of the rotated states |x ± nm ⟩ = 1 2 (|n⟩ ± |m⟩) and |y ± nm ⟩ = 1 2 (|n⟩ ± i|m⟩): Since these involve polulation differences, we can measure them with classical estimators in the same way as we measure populations, where Our mapping for the coherences is therefore which in view of Eq. ( 32) can be written in the neater form Together with the population estimator in Eq. ( 29), this leads to a general estimator for electronic observables that will give the same result in any unitarily rotated basis:

E. Initial conditions
To simulate a non-equilibrium process starting in a pure electronic state |n⟩⟨n|, such as the bright state of a molecule that has just been photoexcited, we need to define an initial electronic distribution, ρ n (c), such that In this way, the initial estimate of the population of state m = n is 1, and that of all other states is zero.There are many choices of ρ n that fulfil this condition.The one we will use in the following is corresponding to a uniform distribution in the region where P n is the largest population (see the coloured regions in Fig. 2).This choice fulfils Eq. ( 36) by the same argument as in Appendix C, regardless of whether the initial state is a diabat or an adiabat.For a twostate problem in which the initial state is specified in the adiabatic basis, Eq. ( 37) is equivalent to the prescription for A(c) that Mannouch and Richardson used in Eq. ( 11).However, our prescription for B(c) is different from theirs, and our formulation avoids their weight function W AB (c).
In practice, one can sample points c from Θ n (c) is as follows.First, sample N pairs (x k , y k ), where x k and y k are Gaussian deviates with zero mean and unit variance.Then set c k = (x k + iy k )/ l (x 2 l + y 2 l ).This has |c| = 1 and, due to the rotational invariance of the multidimensional normal distribution, is a point with uniform probability on the simplex.Finally, check which state has the largest |c k | 2 -if it is the initial state n, then accept the point, and if not, resample a new point in the same way until it is accepted.

IV. APPLICATION TO EXCITON ENERGY TRANSFER
To test the multi-state algorithm introduced in Sec.III, we have applied it to a Frenkel-exciton model of energy transfer in the Fenna-Matthews-Olson complex. 29This has become a standard benchmark for nonadiabatic dynamics and it allows a comparison with a variety of other trajectory-based methods.The problem is challenging for conventional surface hopping due to the presence of 'trivial' crossings with low hopping probabilities, which can require an extremely small time step to converge. 30everal phase-space methods have been reported to perform well for FMO at high temperatures, including spin mapping, 11,31 MMST mapping, [32][33][34][35][36][37] and the symmetric quasi-classical trajectory method of Cotton and Miller. 38owever, these methods struggle to recover the correct long-time equilibrium populations at low temperatures, and they only avoid unstable inverted potentials because all of the potential-energy surfaces in the model are harmonic with the same frequencies.
The standard FMO model is comprised of seven sites, but to facilitate a comparison with a previous surface hopping study we have also considered a reduced model with three sites.In both cases, the sites are coupled independently to identical harmonic baths.The full Hamiltonian and all model parameters are given in Appendix D. The initial condition for our simulations is a pure electronic state in an uncoupled classical Boltzmann bath.The initial electronic state can be in the site or the exciton basis and we shall report results for both below.Each simulation was averaged over 10 5 trajectories for good statistical convergence, although rough results can typically be obtained with 10 3 .We used a time step of 0.25 fs with the trajectory integrator and hopping protocol described in Appendix E. Fully quantum mechanical HEOM benchmark results were computed for comparison using the pyrho open source software. 39 Three-state model Sindhu and Jain have recently used a three-site FMO model to compare a variety of different decoherence corrections to fewest switches surface hopping. 30,40To assess the performance of MASH against FSSH, we shall use the most accurate of the methods they considered in their comparison as a reference.This method, the augmented surface hopping (A-FSSH), employs a parameterfree decoherence correction that has been found to improve upon the original FSSH algorithm for condensedphase problems 17 (for example, it has been found to recover Marcus theory rates in the golden-rule limit 41 ).It should be noted that Sindhu and Jain quantized one of the nuclear modes in their FMO calculation to give a manifold of vibronic states, whereas our MASH calculations treat all nuclear modes classically.Since their quantization improved the agreement of A-FSSH with HEOM rather than harmed it, 30 we still feel this provides a fair comparison.
In Figure 3, we compare MASH to the A-FSSH results from Ref. 30 as well as to HEOM.We find that MASH almost perfectly captures the initial coherent oscillations in the HEOM polulation dynamics, as well as the long-time relaxation.In contrast, A-FSSH leads to overly damped oscillations and too rapid thermalization.These results are interesting not only because they demonstrate that MASH provides a clear improvement over A-FSSH, but because it does so without decoherence corrections.This suggests that overcoherence may not be such a useful way to understand the problems of FSSH as previously thought.The issue may instead be more closely related to the inconsistency between ⟨|c a | 2 ⟩ and the fraction of trajectories on the active surface.This is precisely the problem that MASH was designed to resolve by uniquely defining the active surface in terms of the wavefunction coefficients. 14nother interesting observation is that, even at 77 K where k B T = 53.5 cm −1 is roughly two times smaller than the characteristic phonon energy of ℏω c = 106 cm −1 , MASH gives accurate results despite treating the nuclei classically.This indicates that nuclear quantum effects are less important for this system than the quantization of one of the nuclear modes in Ref. 30 would suggest them to be.

B. Seven-state model
We are not aware of any calculations using conventional surface hopping for the seven-state FMO model, but it is expected that they would suffer from the same difficulties as in the three-state case.In Figure 4, we compare our MASH results for this model directly to those of the quantum HEOM benchmark.Panel (a) shows the population dynamics in the site basis at 300 K, starting from site 1, and panel (b) shows the dynamics in the exciton basis at 300 K, starting from exciton state 1.The exciton basis is defined as the eigenbasis of the system Hamiltonian in Eq. (D2).In both cases, MASH is seen to agree almost perfectly with the fully quantum benchmark.8]42 Panel (c) of Figure 4 shows the same situation as in panel (a), but at 77 K.This low-temperature regime is more challenging for mapping approaches, several of which predict negative populations.However, MASH continues to agree well with HEOM at short times, except perhaps for the small deviation in the transfer between sites 1 and 2, and it gives the correct quantum-mechanical equilibrium populations in the longtime limit.It is worth repeating that, as in the case of the three-state model considered in Figure 3, these longtime populations are obtained correctly without including nuclear quantum effects in the calculation.While nuclear quantum effects have previously been claimed to be important for the functioning of light-harvesting complexes, 43,44 a more recent study has shown that clas- sical nuclei are in fact sufficient to describe realistic models of FMO at room temperature, 45 and our present results suggest that this conclusion extends even to cryogenic temperatures.(Including nuclear quantum effects might perhaps improve on the MASH description of the short-time population transfer between sites 1 and 2 in panel (c) of Figure 4, but this would only be a small correction.)

C. Coherence dynamics
To assess the accuracy of the coherence estimator introduced in Sec.III.D, we have used it to calculate the coherence between sites 1 and 2 of the seven-state FMO model with the initial population on site 1.To facilitate a com- parison with previous spin-mapping simulations, 11 the timescale of the bath was set to τ c = 100 fs in these calculations, as it was in the HEOM calculations of Sarovar et al. 46 .Since these earlier HEOM calculations included a trapping rate that is not present in our FMO model, we recomputed the HEOM results without any trapping, and found that they differ only very slightly (by around 1%) from the results reported in Ref. 46. Figure 5 compares the HEOM coherence with that obtained from MASH.The agreement is not perfect, but it is clearly very good for the imaginary part of the coherence and reasonably good for the real part, especially at short times.Our coherence estimator has therefore been validated for a multi-state problem.There are no previous surface hopping results for this problem, for the same reason as there are none for the population dynamics of the seven-state model in Figure 4.The coherence has, however, been calculated previously using spin mapping, 11 which was found to give results of similar quality to the MASH results in Figure 5.

V. APPLICATION TO TWO-LEVEL SYSTEMS
Since our scheme for measuring electronic observables differs from that of Mannouch and Richardson even for two-level systems, it is important to check the performance of our approach also in this case.For this purpose, we have applied our method to some of the twostate Tully, 1 spin-boson, and pyrazine models considered previously by Mannouch and Richardson. 14All model definitions and parameters are the same as in Ref. 14, where further comparisons to FSSH, Ehrenfest, and spinmapping results can be found.Unless otherwise stated, the results were computed from 10 5 trajectories and integrated with a timestep of 0.25 fs.

A. Tully models
First, in Fig. 6, we consider the problem of wavepacket splitting in Tully's single avoided crossing model. 1,47annouch and Richardson have shown that their scheme reproduces the correct final momentum distribution at both low and high incoming kinetic energies.Our scheme uses a different trajectory weighting that reproduces their results at high energy but gives slightly less accurate peak heights at low energy.The area under each peak is proportional to the fraction of trajectories emerging on each adiabat.In the original MASH, this fraction is ⟨Θ a (t)⟩, where t is a time at which the crossing is complete.In our scheme, the time-dependent populations ⟨Φ a (t)⟩ do give the correct adiabatic state branching ratio in the low energy example.However, since ⟨Φ a (t)⟩ is not the same as the fraction of trajectories emerging on each adiabat, our final momentum distribution is slightly incorrect.Fig. 7 shows the probability of transmission on the upper adiabat for Tully's double avoided crossing model as a function of the initial momentum on the lower adia- bat.Here, our results are similar to those of the original MASH method, and arguably slightly better at high momenta.For low momenta, neither method can reproduce the Stückelberg oscillation, which is due to electronic interference.Note also that our estimator predicts negative transmission probabilities at low momenta where the upper product adiabat is energetically inaccessible, whereas Mannouch and Richardson's transmission probability goes correctly to zero.Despite the issues that these tests have identified, the overall impression we get from the comparisons in Figs. 6  and 7 is that our method does not perform significantly worse for these models than the original MASH method of Mannouch and Richardson. 14Especially when one considers that these non-ergodic, one-dimensional, microcanonical models do not satisfy the assumptions we made when deriving our population estimator, and are not therefore the sort of problems for which our method was designed.The present MASH results in Fig. 6 are certainly more accurate than those of either Ehrenfest dynamics or spin mapping, for example, both of which fail to describe the wavepacket bifurcation. 14

B. Spin-boson model
Next, we consider the spin-boson model, for which where ) is a bath of oscillators with spectral density For this model, we initialized the nuclei from a Wigner distribution to be consistent with Ref. 14.The resulting population dynamics is shown in Fig. 8 for a variety of regimes.The integration time step (in units of ∆ −1 ) was 0.01 for panels (a), (b) and (d), and 0.002 for panel (c).
In scheme is so accurate.The fact that the present version of MASH does not capture the inverted regime correctly is probably the most serious deficiency of the method we have found so far.It is not clear how to solve this problem without abandoning the basis-set independence of our population estimator, but it is conceivable that a quantum-jump procedure 14 might improve our results.

C. Pyrazine
Finally, we consider ultrafast dynamics through the conical intersection in a full-dimentional (24-mode) model of pyrazine.This model includes bilinear couplings which make it challenging for traditional mapping methods.It is also typical of the sort of photochemical problems to which one might expect methods like MASH to be applied.As shown in Fig. 9, our results for this final two-state problem are of comparable quality to those of the original MASH scheme. 14

VI. CONCLUDING REMARKS
In this article, we have shown how to extend the MASH methodology of Mannouch and Richardson 14 to general N -level systems.We have also proposed simplified estimators for populations and coherences that can be applied in any basis and are guaranteed to give the correct long-time equilibrium populations in this basis provided the system is ergodic and the classical approximation to the nuclear motion is valid.Our applications to FMO exciton models have shown that the resulting multi-state MASH method is at least as accurate as any previous semiclassical mapping method, and significantly more accurate than an up-to-date implementation of fewestswitches surface hopping.
Regarding the restriction to classical nuclear motion, we would point out that the majority of interesting nonadiabatic dynamics problems follow a photoexcitation step in which a vast amount of energy (significantly larger than k B T ) has been deposited in the system.The resulting nuclear motion will often be sufficiently fast that the classical nuclear motion approximation is well justified, as it certainly seems to be in all of the FMO calculations we have presented here.However, nuclear quantum effects are expected to be important in some other contexts.For example, nuclear tunnelling is known to have a significant impact on the rates of electron transfer reactions in the Marcus inverted regime. 49It might therefore be interesting to add nuclear quantum effects to the present methodology, perhaps by adapting the ring-polymer molecular dynamics techniques that have already been developed for standard surface hopping. 50egarding how much more successful MASH is for the FMO problem than fewest switches surface hopping (as we have shown in Figure 3), we would say that Mannouch and Richardson's idea of deterministically tying the active surface to the adiabatic state with the largest population was quite inspired.This allowed them to derive the momentum rescaling and momentum reversal stages of their MASH surface hopping algortihm from first principles, to avoid the inconsistency between the stochastically averaged ⟨|c a | 2 ⟩ and the active adiabatic surface V a (q), and thereby to eliminate the need for ad hoc 'decoherence' corrections, in a single stroke.All we have done here is to show that their idea can be adapted to treat an arbitrary number of coupled electronic states in a straightforward and internally consistent way.While this has come at the cost of sacrificing some of the advantages of the original MASH scheme for two-state problems, as we have shown in Figures 6 to 8 and discussed in Sec.V, we feel that the extension to more electronic states is worth this sacrifice because it opens up the possibility of applying the method to a far wider variety of interesting nonadiabatic problems.
where |n⟩⟨n| is a projection onto site n.The frequencies ω j and the coupling coefficients κ j of each bath are determined by the same spectral density which is taken to be of Debye form with a reorganization energy of λ = 35 cm −1 and a characteristic phonon frequency of ω c = 1/τ c .Figures 3 and 4 used τ c = 50 fs (ω c = 106.14cm −1 ), whereas Figure 5 used τ c = 100 fs (ω c = 53.07cm −1 ) for reasons discussed in the text.The Debye bath was discretized into 60 modes per site with the same procedure as described on p. 51 in Ref. 52.

Appendix E: Implementation details
To integrate the MASH equations of motion for a finite time step ∆t, we used a simple velocity-Verlet scheme p ← p − ∇V a (q)∆t/2 (E1b) q ← q + (p/m)∆t (E1c) p ← p − ∇V a (q)∆t/2 (E1d) where a is the index of the active surface.
After each time step, the adiabatic populations were calculated, and if a new state b had reached a higher population than state a, we used 10 bisections to find the crossing point δt < ∆t where P a (δt) = P b (δt).In each bisection iteration, Eq. (E1) was used to propagate the system on state a from the original starting point, but through a time step of δt rather than ∆t.After the hop, the remainder of the original time step, ∆t − δt, was processed in the same way, starting on the new active surface if the hop had been successful.Occasionally, after an unsuccessful hop, we found that the trajectory would attempt to hop again multiple times before reaching the end of the full time step ∆t.There may be better ways to deal with this, but in the present calculations we simply abandoned the trajectory after 30 unsuccessful hopping attempts within any given time step.This happened for less than 1.6% of the trajectories in the worst case.
The hops were implemented by switching the active surface and rescaling the momentum to conserve the total energy, or reversing the momentum and abandoning the hop if the kinetic energy was insufficient to cross the If the argument of the square root is negative, the hop is abandoned and the projection is instead reversed.The resulting rescaling/reversal can be regarded as arising from a classical particle incident on a step barrier, as in the two-level case. 14This way of treating momentum reversal arises naturally from the MASH equations of motion and it differs from the conventional protocols used in FSSH. 53

FIG. 1 .
FIG.1.Bloch sphere representation of a two-level system.Any selection of two opposite fixed points on the sphere (e.g., the black dots) defines a diabatic basis, whereas the adiabatic basis corresponds to the opposite points on an axis that follows (Vx, Vy, Vz) (blue and red dots).In MASH, the nuclei evolve on the adiabatic state that is closest to the instantaneous direction of the Bloch vector, σ = ⟨ψ| σ|ψ⟩.
where h(x) is the Heaviside step function.The weight function W AB (c) is also specified in the adiabatic representation, as W AB (c) = 3 if A and B are both coherences, 2 if one is a coherence and the other a population, and 2|(|c 1 | 2 − |c 2 | 2 )| if they are both populations.Diabatic observables first have to be converted to the adiabatic representation, as in conventional surface hopping, in order to calculate these quantities.

FIG. 3 .
FIG. 3. Comparison of HEOM (dash-dotted), MASH (solid), and A-FSSH (dotted) for population dynamics in a three-site FMO model.Results are reported in the site basis at 77 K, starting in site 1.The A-FSSH results are from Ref. 30.

FIG. 4 .
FIG. 4. Comparison of MASH (solid) and HEOM (dash-dotted) for population dynamics in the seven-site FMO model at (a) 300 K in the site basis, (b) 300 K in the exciton basis, and (c) 77 K in the site basis.The initial state was site 1 for panels (a) and (c) and exciton 1 for panel (b).

FIG. 5 .
FIG.5.Comparison of MASH (solid) and HEOM (dashdotted) results for the dynamics of the dominant coherence of the seven-state FMO model in the site basis at 300 K, starting from site 1.

FIG. 6 .
FIG.6.Final momentum distribution after a single avoided crossing in Tully's first model, for two initial energies.The original MASH results and the model parameters are from Ref. 14.These histograms were generated from 10 6 trajectories.

FIG. 7 .
FIG. 7. Transmission probability from the lower to the upper adiabat in Tully's double avoided crossing model.The original MASH results and the model parameters are from Ref. 14.

1 √
potential step.The time derivative of the adiabatic population difference associated with a hop isṖa − Ṗb = 2 j p j m j a ′ Re c * a ′ d j a ′ a c a − c * a ′ d j a ′ b c b , (E2)from which one can identify the component of the momentum that needs to be rescaled or reversed.In practice, we project the mass-scaled momenta pj = p j / √ m j onto the direction δab = δ ab /|δ ab | of a vector with elementsδ j ab = m j a ′ Re c * a ′ d j a ′ a c a − c * a ′ d j a ′ b c b ,(E3)and rescale the magnitude of this projection from|p old | to |p new | = |p old | 2 + 2(V a − V b ), leaving its orthogonal complement unchanged.