# A multi-state mapping approach to surface hopping

Johan E. Runeson  Department of Chemistry, University of Oxford, Physical and
Theoretical Chemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, UK  David
E. Manolopoulos  Department of Chemistry, University of Oxford, Physical and
Theoretical Chemistry Laboratory, South Parks Road, Oxford, OX1 3QZ, UK

###### Abstract

We describe a multiple electronic state generalisation of the mapping approach
to surface hopping introduced recently by Mannouch and Richardson (J. Chem.
Phys. 158, 104111 (2023)). This generalisation 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 accurate as those of any other semiclassical
method we are aware of.

##  I Introduction

Fewest switches surface hopping (FSSH) is the dominant method for simulating
electronically nonadiabatic dynamics in chemistry. Originally proposed by
Tully in the early 1990’s,Tully (1990) it has become a standard tool for
running ab initio trajectories and is now routinely used to study various
aspects of photochemical dynamics including both energy and charge transfer
processes.Crespo-Otero and Barbatti (2018); Menger  _et al._ (2018); Peng  _et
al._ (2022); Toldo  _et al._ (2023) 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,Subotnik, Ouyang, and Landry (2013) 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,Wang  _et al._ (2020) 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.Parandekar and Tully (2006) Semiclassical mapping methods represent
the multiple electronic states as harmonic oscillatorsMeyer and Miller (1979);
Stock and Thoss (1997) or generalized spins,Runeson and Richardson (2020) and
use either full phase space sampling or symmetrical quasi-classical
binningCotton and Miller (2013); Miller and Cotton (2016) 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 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).Mannouch and Richardson (2023) As 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 procedureKapral (2015)). And 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.Mannouch and Richardson (2023) It
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.Mannouch and Richardson (2023) It
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 generalization 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 generalization 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,Mannouch and Richardson (2023) and in Sec.
V we conclude this paper.

##  II Nonadiabatic dynamics

Consider a general nonadiabatic system defined by the Hamiltonian

$\hat{H}(p,q)=\sum_{j=1}^{f}\frac{p_{j}^{2}}{2m_{j}}+\hat{V}(q),$ (1)

where $q=\\{q_{j}\\}$ and $p=\\{p_{j}\\}$ are the coordinates and momenta of
$f$ nuclear degrees of freedom with masses $\\{m_{j}\\}$, and $\hat{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\rangle\\}$, i.e.,

$\hat{V}(q)=\sum_{nm}V_{nm}(q)|n\rangle\langle m|.$ (2)

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 $\hat{V}(q)$), where

$\hat{V}(q)=\sum_{a}V_{a}(q)|a(q)\rangle\langle a(q)|.$ (3)

The two bases are related through their transformation matrix elements
$U_{na}(q)=\langle n|a(q)\rangle$. 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 long-standing problem
in semiclassical dynamics. To see why, consider a particle at configuration
$q$ with momentum $p$ and with electronic state
$|\psi\rangle=\sum_{n}c_{n}|n\rangle$. 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

$\dot{c}_{n}=-\mathrm{i}\sum_{m}V_{nm}(q)c_{m},$ (4)

or in the adiabatic basis

$\dot{c}_{a}=-\mathrm{i}V_{a}(q)c_{n}-\sum_{j}\frac{p_{j}}{m_{j}}\sum_{b}d_{ab}^{j}(q)c_{b},$
(5)

where $d^{j}_{ab}(q)=\langle a(q)|\nabla_{j}|b(q)\rangle$ is an element of the
nonadiabatic coupling vector. (Throughout this paper we use units where
$\hbar=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’ $\hat{F}_{j}(q)=-\nabla_{j}\hat{V}(q)$,
but this is not yet useful to run classical dynamics. What one would like is
equations of motion of the form

$\displaystyle\dot{q}_{j}$ $\displaystyle=p_{j}/m_{j}$ (6a)
$\displaystyle\dot{p}_{j}$ $\displaystyle=F_{j}(q),$ (6b)

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.

###  II.1 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)=-\langle a(q)|\nabla_{j}\hat{V}|a(q)\rangle.$ (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

$P_{a\to b}=2\Delta
t\operatorname{Re}\left(\frac{c_{b}}{c_{a}}\right)\sum_{j}\frac{p_{j}}{m_{j}}d_{ab}^{j}.$
(8)

Tully chose this probability such that the fraction of trajectories evolving
on surface $a$ would approximate the average population
$\langle|c_{a}|^{2}\rangle$ with a minimal number of switches.Tully (1990) If
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
Martens (2016)) 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
$\langle|c_{a}|^{2}\rangle$. This inconsistency is at the root of many of the
issues present in surface hopping.Subotnik  _et al._ (2016); Carof, Giannini,
and Blumberger (2017) Traditionally, 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 long-time equilibrium
in condensed-phase systems, although it does often provide a better
approximation than Ehrenfest dynamics.Parandekar and Tully (2005); Schmidt,
Parandekar, and Tully (2008)

###  II.2 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,

$F_{j}(q)=-\langle\psi|\nabla_{j}\hat{V}(q)|\psi\rangle=-\sum_{nm}\nabla_{j}V_{nm}(q)c_{n}^{*}c_{m}.$
(9)

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. Meyer and Miller (1979); Stock and
Thoss (1997) Taking 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.Stock and Thoss (2005); Miller and
Cotton (2016)

Harmonic oscillators are not the only way to map the electronic coefficients
onto classical variables. For two-level systems, another choice would be to
use the well-known isomorphism to a spin-1/2 system. This approach was
recently used to develop a ‘spin mapping’ analogous to the MMST
mapping,Runeson and Richardson (2019) 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.Runeson and Richardson (2020) This 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 $\beta=1/k_{\rm
B}T$.Runeson  _et al._ (2022a)

However, just like the MMST mapping, spin mapping can predict (unphysical)
negative populations for 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.Amati,
Runeson, and Richardson (2023) This 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.

###  II.3 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.

Figure 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 $(V_{x},V_{y},V_{z})$ (blue and red dots). In MASH,
the nuclei evolve on the adiabatic state that is closest to the instantaneous
direction of the Bloch vector,
$\bm{\sigma}=\langle\psi|\hat{\bm{\sigma}}|\psi\rangle$.

Mannouch and Richardson have recently shown that, at least for two-level
systems, this may indeed be possible.Mannouch and Richardson (2023) Their
‘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 $a$ separately from the electronic wavefunction, they set
it to be the 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. Moreover, the inconsistency between the active surface $a$ and
the stochastically averaged probability $\langle|c_{a}|^{2}\rangle$ in
traditional surface hopping is eliminated entirely.

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 $\bm{V}=(V_{x},V_{y},V_{z})$, where
$V_{\alpha}=\mathrm{Tr}[\hat{V}\hat{\sigma}_{\alpha}]$ and the
$\\{\hat{\sigma}_{\alpha}\\}$ are the Pauli spin operators. In MASH, the
instantaneous active surface is set to 1 when the Bloch vector
$\bm{\sigma}=\langle\psi|\hat{\bm{\sigma}}|\psi\rangle$ 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 solving the inconsistency problem, the MASH approach has
several other appealing features, a more detailed discussion of which can be
found in Ref. Mannouch and Richardson, 2023:

  1. (i)

it satisfies detailed balance in the sense that it gives the correct long-time
populations of the adiabatic electronic states (assuming ergodicity);

  2. (ii)

the momentum rescaling and momentum reversal arise naturally in the
deterministic surface hopping algorithm without any ambiguity;

  3. (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.Runeson and Richardson (2021) They also found MASH
to be the most accurate classical-trajectory method considered so far in
describing ultrafast internal conversion in pyrazine.Mannouch and Richardson
(2023)

However, their method is currently limited to two-level systems and it is not
obvious how to generalize it to more electronic states, partly because
Mannouch and Richardson specifically constructed their observables for the
case of two. Explicitly, for electronic observables $\hat{A}$ and $\hat{B}$,
they computed the time-correlation function 111In their paper, Mannouch and
Richardson used spin vectors rather than $c$ coordinates, but these contain
the same information and are simply a different a notation for the same
quantity.

$C_{AB}(t)=\langle\mathcal{W}_{AB}(c_{0})A(c_{0})B(c_{t})\rangle,$ (10)

where $\langle\cdot\rangle$ 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,

$\displaystyle|1\rangle\langle 1|$ $\displaystyle\mapsto
h(|c_{1}|^{2}-|c_{2}|^{2})$ (11a) $\displaystyle|2\rangle\langle 2|$
$\displaystyle\mapsto h(|c_{2}|^{2}-|c_{1}|^{2})$ (11b)
$\displaystyle|a\rangle\langle b|$ $\displaystyle\mapsto c_{a}^{*}c_{b}\quad\
\,(a\neq b),$ (11c)

where $h(x)$ is the Heaviside step function. The weight function
$\mathcal{W}_{AB}(c)$ is also specified in the adiabatic representation, as
$\mathcal{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.

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 populations 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 extension 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 extension to 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.

###  III.1 Dynamics

The dynamics extends quite naturally from the two-level 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

$\displaystyle P_{n}$ $\displaystyle=|c_{n}|^{2}$ (12a)
$\displaystyle\phi_{n}$ $\displaystyle=\arg c_{n}.$ (12b)

If we also define the classical state projectors

$\Theta_{n}=\begin{cases}1&\text{if }P_{n}\geq P_{m}~{}\forall\,m\neq n\\\
0&\text{otherwise},\end{cases}$ (13)

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 adiabatic basis
state $a$ with a non-zero classical state projector $\Theta_{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,
$\sum_{n=1}^{N}P_{n}=\sum_{n=1}^{N}|c_{n}|^{2}=1$, and they will always be
non-negative, $P_{n}\geq 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},\dots,P_{N})$ is a point on
the _simplex_

$\\{P\in\mathbb{R}^{N}|\sum_{n}P_{n}=1,P_{n}\geq 0~{}\forall\,n\\}$ (14)

which is a geometrical object that can be thought of as a high-dimensional
triangle with vertices on the unit vectors.

Figure 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.

When 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 the nonadiabatic coupling vector. If the kinetic energy is
insufficient to hop, we reverse the momentum along this vector. So far,
everything is identical to the two-level case.Mannouch and Richardson (2023)
In 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-mechanical equilibrium populations in the long-
time limit, we shall begin by discussing what happens when the dynamics has
reached equilibrium.

###  III.2 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 long-time limit reduce to the equilibrium
distribution

$\rho(q,p,c)\propto\mathrm{e}^{-\beta E(p,q,c)},$ (15)

where the energy function is

$E(p,q,c)=T(p)+\sum_{a}V_{a}(q)\Theta_{a}$ (16)

with $T(p)=\sum_{j}p_{j}^{2}/2m_{j}$ and $\Theta_{a}\equiv\Theta_{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.

For a coupled electronic–nuclear system, the quantum-mechanical partition
function with classical nuclei is

$Z=\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\,\mathrm{Tr}[\mathrm{e}^{-\beta\hat{H}(p,q)}].$
(17)

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$,

$Z_{\rm
MASH}=\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\int_{|c|=1}\frac{\mathrm{d}c}{\mathcal{N}}\,\mathrm{e}^{-\beta
E(p,q,c)},$ (18)

where $\mathcal{N}=2\pi^{N}/N!$ normalizes the electronic integral such that
$\int_{|c|=1}{\mathrm{d}c}/{\mathcal{N}}=N$ (see Appendix A). Because each
$\Theta_{a}$ is an idempotent projector onto a different region of $c$ space,
and it has a unit integral
$\int_{|c|=1}{\mathrm{d}c}/{\mathcal{N}}\,\Theta_{a}=1$, it follows that

$\int_{|c|=1}\frac{\mathrm{d}c}{\mathcal{N}}\,\mathrm{e}^{-\beta\sum_{a}V_{a}(q)\Theta_{a}}=\sum_{a}\mathrm{e}^{-\beta
V_{a}(q)},$ (19)

which shows that MASH is consistent with the mixed quantum-classical partition
function ($Z_{\rm MASH}=Z$). By the same argument, the thermal equilibrium
population of each adiabatic quantum state,

$\langle{P}_{a}\rangle=\frac{1}{Z}\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\,\mathrm{Tr}[\mathrm{e}^{-\beta\hat{H}(p,q)}|a(q)\rangle\langle
a(q)|],$ (20)

will be equal to the corresponding MASH expression

$\langle\Theta_{a}\rangle=\frac{1}{Z}\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\int_{|c|=1}\frac{\mathrm{d}c}{\mathcal{N}}\,\mathrm{e}^{-\beta
E(p,q,c)}\Theta_{a},$ (21)

provided $\Theta_{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 $\langle\Theta_{n}\rangle$
does not equal $\langle P_{n}\rangle$ in other bases. To see this, we can
write the latter explicitly as

$\displaystyle\langle P_{n}\rangle$
$\displaystyle=\frac{1}{Z}\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\,\mathrm{Tr}[\mathrm{e}^{-\beta\hat{H}(p,q)}|n\rangle\langle
n|]$
$\displaystyle=\frac{1}{Z}\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\,\sum_{a}\mathrm{e}^{-\beta[T(p)+V_{a}(q)]}|\langle
a(q)|n\rangle|^{2},$ (22)

whereas the former is

$\displaystyle\langle\Theta_{n}\rangle$
$\displaystyle=\frac{1}{Z}\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\int_{|c|=1}\frac{\mathrm{d}c}{\mathcal{N}}\,\mathrm{e}^{-\beta
E(p,q,c)}\Theta_{n}$
$\displaystyle=\frac{1}{Z}\int\frac{\mathrm{d}p\mathrm{d}q}{(2\pi)^{f}}\sum_{a}\mathrm{e}^{-\beta[T(p)+V_{a}(q)]}\int_{|c|=1}\frac{\mathrm{d}c}{\mathcal{N}}\,\Theta_{a}\Theta_{n}.$
(23)

This is in general _not_ equal to $\langle P_{n}\rangle$ because

$\int\frac{\mathrm{d}c}{\mathcal{N}}\,\Theta_{a}\Theta_{n}\neq|\langle
a(q)|n\rangle|^{2}.$ (24)

In other words, the state projection defined by Eq. (13) does not transform
correctly as a population estimator under unitary basis rotations.

###  III.3 Populations

A population estimator $\Phi_{n}$ can in fact be constructed such that
$\langle\Phi_{n}\rangle=\langle P_{n}\rangle$ 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 a population estimator of the form

$|n\rangle\langle n|\mapsto\Phi_{n}=\alpha_{N}\,P_{n}+\beta_{N},$ (25)

in which $\alpha_{N}$ and $\beta_{N}$ are scalar parameters that are to yet be
determined. Our goal is to choose these parameters such that

$\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}\Phi_{n}=|\langle
a(q)|n\rangle|^{2}$ (26)

holds for all adiabatic states $a$ and diabatic states $n$. This will suffice
to ensure – by virtue of the argument in Eqs. (III.2) and (III.2) – that the
expectation value of $\Phi_{n}$ will be correct at equilibrium. Appendix B
shows that Eq. (26) will be satisfied if the simpler condition

$\int_{|c|=1}\frac{\mathrm{d}c}{\mathcal{N}}\,\Theta_{a}\Phi_{b}=\delta_{ab}$
(27)

holds for all adiabatic states $a$ and $b$. The constraints that this imposes
on $\alpha_{N}$ and $\beta_{N}$ are then investigated in Appendix C, which
obtains the unique solution

$\alpha_{N}={N-1\over H_{N}-1}\quad\hbox{and}\quad\beta_{N}={1-\alpha_{N}\over
N},$ (28)

where $H_{N}=\sum_{n=1}^{N}1/n$. Hence we arrive at the estimator

$\Phi_{n}={1\over N}+\alpha_{N}\left(P_{n}-{1\over N}\right),$ (29)

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).

###  III.4 Coherences

In contrast to Ref. Mannouch and Richardson, 2023, as well as most other
literature on surface hopping, we shall rewrite 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
we have shown that the population estimator in Eq. (29) transforms correctly
under unitary basis rotations.

The coherences between states $|n\rangle$ and $|m\rangle$ are linear
combinations of the Hermitian operators
$\hat{\sigma}^{x}_{nm}=|n\rangle\langle m|+|m\rangle\langle n|$ and
$\hat{\sigma}^{y}_{nm}=-\mathrm{i}\left(|n\rangle\langle m|-|m\rangle\langle
n|\right)$,

$\displaystyle|n\rangle\langle m|$ $\displaystyle={1\over
2}\left(\hat{\sigma}^{x}_{nm}+\mathrm{i}\,\hat{\sigma}^{y}_{nm}\right)$ (30a)
$\displaystyle|m\rangle\langle n|$ $\displaystyle={1\over
2}\left(\hat{\sigma}^{x}_{nm}-\mathrm{i}\,\hat{\sigma}^{y}_{nm}\right),$ (30b)

and we can rewrite $\hat{\sigma}^{x}_{nm}$ and $\hat{\sigma}^{y}_{nm}$ as
differences between the population operators of the rotated states
$|x^{\pm}_{nm}\rangle=\sqrt{1\over 2}\left(|n\rangle\pm|m\rangle\right)$ and
$|y^{\pm}_{nm}\rangle=\sqrt{1\over
2}\left(|n\rangle\pm\mathrm{i}|m\rangle\right)$:

$\displaystyle\phantom{1\over 2}\hat{\sigma}^{x}_{nm}$
$\displaystyle=|x^{+}_{nm}\rangle\langle x^{+}_{nm}|-|x^{-}_{nm}\rangle\langle
x^{-}_{nm}|$ (31a) $\displaystyle\phantom{1\over 2}\hat{\sigma}^{y}_{nm}$
$\displaystyle=|\,y^{+}_{nm}\rangle\langle
y^{+}_{nm}|-|\,y^{-}_{nm}\rangle\langle y^{-}_{nm}|.$ (31b)

Since these involve polulation differences, we can measure them with classical
estimators in the same way as we measure populations,

$\displaystyle\phantom{1\over
2}\hat{\sigma}^{x}_{nm}\mapsto\Sigma^{x}_{nm}=\alpha_{N}\left(P_{x^{+}_{nm}}-P_{x^{-}_{nm}}\right)\phantom{,}$
(32a) $\displaystyle\phantom{1\over
2}\hat{\sigma}^{y}_{nm}\mapsto\Sigma^{y}_{nm}=\alpha_{N}\left(P_{y^{+}_{nm}}-P_{y^{-}_{nm}}\right),$
(32b)

where $P_{x^{\pm}_{nm}}=|c_{x^{\pm}_{nm}}|^{2}$ and
$P_{y^{\pm}_{nm}}=|c_{y^{\pm}_{nm}}|^{2}$ with $c_{x^{\pm}_{nm}}=\sqrt{1\over
2}(c_{n}\pm c_{m})$ and $c_{y^{\pm}_{nm}}=\sqrt{1\over
2}(c_{n}\pm\mathrm{i}\,c_{m})$. Our mapping for the coherences is therefore

$\displaystyle|n\rangle\langle m|$ $\displaystyle\mapsto{1\over
2}\left(\Sigma^{x}_{nm}+\mathrm{i}\,\Sigma^{y}_{nm}\right)$ (33a)
$\displaystyle|m\rangle\langle n|$ $\displaystyle\mapsto{1\over
2}\left(\Sigma^{x}_{nm}-\mathrm{i}\,\Sigma^{y}_{nm}\right).$ (33b)

We would stress again that this automatically treats the coherences on an
equal footing to the populations: coherences are no more “quantum mechanical”
(and therefore no harder to treat classically) within this framework than
populations. All that matters is that the population estimator that is used
can be applied equally well in any unitarily rotated basis, as ours can.

###  III.5 Initial conditions

To simulate a non-equilibrium process starting in a pure electronic state
$|n\rangle\langle n|$, such as the bright state of a molecule that has just
been photoexcited, it is natural to pick an initial value for the electronic
state vector $c$ such that our population estimators satisfy
$\Phi_{m}(c)=\delta_{mn}$. In view of Eq. (29), this implies that

$|c_{m}|^{2}=P_{m}={1\over\alpha_{N}}\delta_{nm}+{\alpha_{N}-1\over
N\alpha_{N}}.$ (34)

Since the phases $\phi_{m}=\arg c_{m}$ are not determined by this constraint,
we sample them uniformly. The resulting initial conditions are reminiscent of
the so-called ‘focused’ initial conditions used in some mapping
approaches,Bonella and Coker (2003); Müller and Stock (1999) and one may think
of the sampled points as circles on the Bloch sphere.Runeson and Richardson
(2019) The radii of these circles depend on $N$, as they do in spin mapping,
but they have different values. In the two-level case, the radius would
actually be equivalent to using focused initial conditions for the MMST
mapping variables.Mannouch and Richardson (2023)

##  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.Ishizaki and Fleming (2009) This 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.Sindhu and Jain (2022) Several phase-space methods have been reported
to perform well for FMO at high temperatures, including spin mapping,Runeson
and Richardson (2020); Mannouch and Richardson (2020) MMST mapping,Tao and
Miller (2010); Kelly and Rhee (2011); Kim and Rhee (2014); Kim, Lee, and Rhee
(2014); Saller, Kelly, and Richardson (2020); Lee, Huo, and Coker (2016) and
the symmetric quasi-classical trajectory method of Cotton and Miller.Cotton
and Miller (2016) However, 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. We ran 105 trajectories in each of our
simulations, with a time step of 20 a.u. (roughly
$0.5\text{\,}\mathrm{f}\mathrm{s}$). The time step was shortened around
hopping events to improve energy conservation as described in Appendix E.
Fully quantum mechanical HEOM benchmark results were computed for comparison
using the pyrho open source software.Berkelbach (2020)

###  IV.1 Three-state model

Figure 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\text{\,}\mathrm{K}$, starting in site 1. The A-FSSH results
are from Ref. Sindhu and Jain, 2022.

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

Sindhu and Jain have recently used a three-site FMO model to compare a variety
of different decoherence corrections to fewest switches surface hopping.Sindhu
and Jain (2022, 2023) To 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
parameter-free decoherence correction that has been found to improve upon the
original FSSH algorithm for condensed-phase problemsSubotnik  _et al._ (2016)
(for example, it has been found to recover Marcus theory rates in the golden-
rule limitJain and Subotnik (2015)). 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,Sindhu and Jain (2022) we still feel this provides
a fair comparison.

In Figure 3, we compare MASH to the A-FSSH results from Ref. Sindhu and Jain,
2022 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 $\langle|c_{a}|^{2}\rangle$ and the fraction of
trajectories on the active surface. This is precisely the problem that MASH
was designed to solve by uniquely defining the active surface in terms of the
wavefunction coefficients.Mannouch and Richardson (2023)

Another interesting observation is that, even at $77\text{\,}\mathrm{K}$ where
$k_{\mathrm{B}}T=$53.5\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ is roughly two
times smaller than the characteristic phonon energy of $\hbar\omega_{\rm
c}=$106\text{\,}\mathrm{c}\mathrm{m}^{-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. Sindhu and Jain, 2022 would suggest them to be. It is
also easier to deal with ‘trivial’ crossings within MASH, because the
deterministic algorithm allows us to adjust the timestep dynamically (see
Appendix E). In other words, MASH is not only more accurate but also
numerically more tractable than conventional surface hopping.

###  IV.2 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\text{\,}\mathrm{K}$, starting from site 1, and panel (b) shows the
dynamics in the exciton basis at $300\text{\,}\mathrm{K}$, starting from
exciton state 1. The exciton basis is defined as the eigenbasis of the system
Hamiltonian in Eq. (50). In both cases, MASH is seen to agree almost perfectly
with the fully quantum benchmark. The agreement is at least as good as or
better than that obtained with previously reported mapping approaches at the
same temperature.Cotton and Miller (2016); Lee, Huo, and Coker (2016); Liu and
Hanna (2018); Runeson and Richardson (2020); Mannouch and Richardson (2020);
Saller, Kelly, and Richardson (2020)

Panel (c) of Figure 4 shows the same situation as in panel (a), but at
$77\text{\,}\mathrm{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 long-time limit. It
is worth repeating that, as in the case of the three-state model considered in
Figure 3, these long-time 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,Cao  _et al._ (2020); Kundu, Dani, and Makri (2022) a
more recent study has shown that classical nuclei are in fact sufficient to
describe realistic models of FMO at room temperature,Runeson  _et al._ (2022b)
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.)

Figure 5: Comparison of MASH (solid) and HEOM (dash-dotted) results for the
dynamics of the dominant coherence of the seven-state FMO model in the site
basis at $300\text{\,}\mathrm{K}$, starting from site 1.

###  IV.3 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
comparison with previous spin-mapping simulations,Runeson and Richardson
(2020) the timescale of the bath was set to $\tau_{\rm
c}=$100\text{\,}\mathrm{f}\mathrm{s}$$ in these calculations, as it was in the
HEOM calculations of Sarovar et al.Sarovar  _et al._ (2010). 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. Sarovar  _et al._ , 2010.

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,Runeson and Richardson (2020) which was found to
give results of similar quality to the MASH results in Figure 5.

##  V Concluding remarks

In this article, we have shown how to extend the MASH methodology of Mannouch
and RichardsonMannouch and Richardson (2023) 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 fewest-
switches 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_{\rm 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.Lawrence  _et al._ (2019) It
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.Shushkov, Li, and Tully (2012)

Regarding how much more successful MASH is for the FMO problem than fewest
switches surface hopping (as we have shown in Figure 3), we would simply say
that Mannouch and Richardson’s idea of deterministically tying the active
surface to the adiabatic state with the largest population was truly inspired.
This allowed them to derive the momentum rescaling and momentum reversal
stages of their MASH surface hopping algortihm from first principles, to
remove the inconsistency between the stochastically averaged
$\langle|c_{a}|^{2}\rangle$ and the active adiabatic surface $V_{a}(q)$, and
thereby to eliminate the need for any ad hoc ‘decoherence’ corrections, in a
single stroke. All we have done here is to show that MASH can be extended to
treat an arbitrary number of coupled electronic states in an equally
straightforward and internally consistent way.

## Acknowledgements

We are grateful to Oliver Riordan for providing us with the proof that
$\pi_{N}=H_{N}/N$ in Appendix C. Johan Runeson is supported by a mobility
fellowship from the Swiss National Science Foundation.

## Author declarations

### Conflict of interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available within the
article.

##  Appendix A Integrals over quantum states

Our integrals over electronic states are integrals over the real and imaginary
parts of each $c_{n}=x_{n}+\mathrm{i}y_{n}$ subject to the constraint that
$|c|=1$. They can be written in unnormalised form as

$\int_{|c|=1}\mathrm{d}c\,f=\int\prod_{n}(\mathrm{d}x_{n}\mathrm{d}y_{n})\,\delta\left(r-1\right)\,f,$
(35)

where $r=[\,\sum_{n}(x_{n}^{2}+y_{n}^{2})\,]^{1/2}$ and $f=f(c)$. The
transformation $x_{n}=P_{n}^{1/2}\cos\phi_{n}$ and
$y_{n}=P_{n}^{1/2}\sin\phi_{n}$ converts this to an integral over populations
and phases,

$\int_{|c|=1}\mathrm{d}c\,f=2\int_{\sum_{n}P_{n}=1}\prod_{n}\left(\frac{1}{2}\mathrm{d}P_{n}\mathrm{d}\phi_{n}\right)\,f,$
(36)

where the first factor of 2 comes from squaring the constraint,
$\int\mathrm{d}r\,\delta(r-1)=2\int\mathrm{d}r\,\delta(r^{2}-1)$, and the
factors of 1/2 come from the Jacobian of the transformation. When the
integrand is purely a function of the populations, $f=f(P)$, the phases can be
integrated out to leave

$\int_{|c|=1}\mathrm{d}c\,f=2\pi^{N}\int_{\sum_{n}P_{n}=1}\mathrm{d}P\,f,$
(37)

and when $f=1$ the population integral gives the volume $1/(N-1)!$ of the
simplex in Eq. (14):

$\int_{|c|=1}\mathrm{d}c={2\pi^{N}\over(N-1)!}.$ (38)

To ensure the correct trace of the identity operator on the electronic space,
we choose to normalise the integrals such that

$\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}=N!\int_{\sum_{n}P_{n}=1}\mathrm{d}P=N,$
(39)

which gives $\mathcal{N}=2\pi^{N}/N!$.

##  Appendix B An alternative constraint

The goal here is to show that Eq. (26) will be satisfied for all $a$ and $n$
provided Eq. (27) is satisfied for all $a$ and $b$. From the definition of
$\Phi_{n}$ as $\alpha_{N}P_{n}+\beta_{N}$, and the fact that
$\int_{|c|=1}{\mathrm{d}c/\mathcal{N}}\,\Theta_{a}=1$, we can start by writing

$\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}\Phi_{n}=\beta_{N}+\alpha_{N}\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}P_{n}.$
(40)

Now transforming $P_{n}$ into the adiabatic representation gives

$P_{n}=|c_{n}|^{2}=\sum_{bb^{\prime}}c_{b}^{*}\langle b(q)|n\rangle\langle
n|b^{\prime}(q)\rangle c_{b^{\prime}}.$ (41)

When this is substituted into the integral on the right hand side of Eq. (B1),
only the diagonal terms (with $b^{\prime}=b$ and
$c_{b}^{*}c_{b^{\prime}}=|c_{b}|^{2}=P_{b}$) survive, because the phase
factors in the off-diagonal terms integrate to zero. Hence

$\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}P_{n}=\sum_{b}|\langle
b(q)|n\rangle|^{2}\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}P_{b}.$
(42)

Substituting this back into Eq. (B1) and rearranging, using the definition of
$\Phi_{b}$ as $\alpha_{N}P_{b}+\beta_{N}$ and the fact that $\sum_{b}|\langle
b(q)|n\rangle|^{2}=1$, gives

$\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}\Phi_{n}=\sum_{b}|\langle
b(q)|n\rangle|^{2}\int_{|c|=1}{\mathrm{d}c\over\mathcal{N}}\,\Theta_{a}\Phi_{b}.$
(43)

If Eq. (27) is satisfied, this reduces to $|\langle a(q)|n\rangle|^{2}$, which
completes what we set out to show.

##  Appendix C Derivation of the population estimator

We now turn to the problem of finding values for $\alpha_{N}$ and $\beta_{N}$
such that the constraints in Eq. (27) are satisfied for all $a$ and $b$. This
is possible because we can rewrite all $N^{2}$ constraints in terms of the
single parameter

$\int_{|c|=1}{{\rm d}c\over\mathcal{N}}\,\Theta_{a}P_{a}=\pi_{N},$ (44)

which is the same for all $a$ by symmetry. Indeed, since
$\sum_{b=1}^{N}P_{b}=1$ we have that

$\int_{|c|=1}{{\rm d}c\over\mathcal{N}}\,\Theta_{a}P_{b}={N\pi_{N}-1\over
N-1}\,\delta_{ab}-{\pi_{N}-1\over N-1},$ (45)

and since $\int_{|c|=1}{\mathrm{d}c/\mathcal{N}}\,\Theta_{a}=1$ we have that

$\displaystyle\int_{|c|=1}{{\rm d}c\over\mathcal{N}}\,\Theta_{a}\Phi_{b}$
$\displaystyle=\alpha_{N}\left({N\pi_{N}-1\over N-1}\right)\delta_{ab}$
$\displaystyle-\alpha_{N}\left({\pi_{N}-1\over N-1}\right)+\beta_{N},$ (46)

which reduces to $\delta_{ab}$ when we set

$\alpha_{N}={N-1\over N\pi_{N}-1}\quad\hbox{and}\quad\beta_{N}={\pi_{N}-1\over
N\pi_{N}-1}.$ (47)

It remains to calculate $\pi_{N}$, which can be interpreted in view of its
definition in Eq. (C1) as the expectation value of the largest coordinate on
the simplex in Eq. (14). This is the same (by symmetry) as the average value
of $P_{1}$ in the region where $P_{1}\geq P_{2}\geq\ldots\geq P_{N}\geq 0$ and
$\sum_{n=1}^{N}P_{n}=1$, which can found as follows. Let $P_{N+1}=0$ and
consider the average value of $\delta P_{n}=P_{n}-P_{n+1}$ in the region where
$\delta P_{n}\geq 0$ and $\sum_{n=1}^{N}n\,\delta P_{n}=1$. This is just a
shrunken simplex with vertices at $(0,\ldots,0,1/n,0,\ldots,0)$, where $1/n$
is in the $n$-th position. So we can stretch each coordinate $\delta P_{n}$ by
a factor of $n$ to recover a standard simplex. The average value of any
coordinate in a standard simplex is $1/N$, so the average value of $\delta
P_{n}$ in the shrunken simplex is $1/(nN)$, and since
$P_{1}=\sum_{n=1}^{N}\delta P_{n}$ this gives

$\pi_{N}={H_{N}\over N},$ (48)

where $H_{N}=\sum_{n=1}^{N}1/n$ is the $N$-th harmonic number.

##  Appendix D FMO model

The FMO model is a standard Frenkel-exciton Hamiltonian in the single-
excitation manifold,

$\hat{H}=\hat{H}_{\rm s}+\hat{H}_{\rm b}+\hat{H}_{\rm sb}.$ (49)

The seven-site system Hamiltonian $\hat{H}_{\rm s}$ is given in the site basis
in units of cm-1 as the $7\times 7$ matrixAdolphs and Renger (2006)

$\begin{pmatrix}200&-87.7&5.5&-5.9&6.7&-13.7&-9.9\\\
-87.7&320&30.8&8.2&0.7&11.8&4.3\\\ 5.5&30.8&0&-53.5&-2.2&-9.6&6.0\\\
-5.9&8.2&-53.5&110&-70.7&-17.0&-63.3\\\ 6.7&0.7&-2.2&-70.7&270&81.1&-1.3\\\
-13.7&11.8&-9.6&-17.0&81.1&420&39.7\\\
-9.9&4.3&6.0&-63.3&-1.3&39.7&230\end{pmatrix},$ (50)

where we have subtracted the lowest bacteriochlorophyll excitation energy of
$12\,210\text{\,}\mathrm{c}\mathrm{m}^{-1}$ from the diagonal entries to save
paper. The three-site model uses the upper left $3\times 3$ submatrix.

The bath Hamiltonian is

$\hat{H}_{\text{b}}=\sum_{n=1}^{N}\sum_{j=1}^{f}\left(\frac{1}{2}p_{j,n}^{2}+\frac{1}{2}\omega_{j}^{2}q_{j,n}^{2}\right),$
(51)

and the system–bath coupling is

$\hat{H}_{\rm
sb}=\sum_{n=1}^{N}\sum_{j=1}^{f}\kappa_{j}q_{j,n}|n\rangle\langle n|,$ (52)

where $|n\rangle\langle n|$ is a projection onto site $n$. The frequencies
$\omega_{j}$ and the coupling coefficients $\kappa_{j}$ of each bath are
determined by the same spectral density

$J(\omega)=\frac{\pi}{2}\sum_{j}\frac{\kappa_{j}^{2}}{\omega_{j}}\delta(\omega-\omega_{j}),$
(53)

which is taken to be of Debye form

$J(\omega)=2\lambda\frac{\omega\omega_{\rm c}}{\omega^{2}+\omega^{2}_{\rm c}}$
(54)

with a reorganization energy of
$\lambda=$35\text{\,}\mathrm{c}\mathrm{m}^{-1}$$ and a characteristic phonon
frequency of $\omega_{\rm c}=1/\tau_{\rm c}$. Figures 3 and 4 used $\tau_{\rm
c}=$50\text{\,}\mathrm{f}\mathrm{s}$$ ($\omega_{\rm
c}=$106.14\text{\,}\mathrm{c}\mathrm{m}^{-1}$$), whereas Figure 5 used
$\tau_{\rm c}=$100\text{\,}\mathrm{f}\mathrm{s}$$ ($\omega_{\rm
c}=$53.07\text{\,}\mathrm{c}\mathrm{m}^{-1}$)$ for reasons discussed in the
text. The Debye bath was discretized into 60 modes per site with the same
procedure as described in Ref. Hele, 2011.

##  Appendix E Implementation details

To integrate the MASH equations of motion for a finite time step $\Delta t$,
we used a simple velocity-Verlet scheme

$\displaystyle c_{n}$
$\displaystyle\leftarrow\sum_{m}[\mathrm{e}^{-\mathrm{i}\hat{V}(q)\Delta
t/2}]_{nm}c_{m}$ (55a) $\displaystyle q$ $\displaystyle\leftarrow q+\Delta
t\,p/2m$ (55b) $\displaystyle p$ $\displaystyle\leftarrow p-\Delta t\,\nabla
V_{a}(q)$ (55c) $\displaystyle q$ $\displaystyle\leftarrow q+\Delta t\,p/2m$
(55d) $\displaystyle c_{n}$
$\displaystyle\leftarrow\sum_{m}[\mathrm{e}^{-\mathrm{i}\hat{V}(q)\Delta
t/2}]_{nm}c_{m}$ (55e)

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
higher a population than state $a$, the active surface was switched and the
momentum was rescaled to conserve the total energy (or reversed and the hop
abandoned if the kinetic energy was insufficient to cross the potential step).
Since the adiabatic population difference associated with the hop has the time
derivative

$\dot{P}_{a}-\dot{P}_{b}=2\sum_{j}\frac{p_{j}}{m_{j}}\sum_{a^{\prime}}\operatorname{Re}\left[c_{a^{\prime}}^{*}d^{j}_{a^{\prime}a}c_{a}-c_{a^{\prime}}^{*}d^{j}_{a^{\prime}b}c_{b}\right],$
(56)

the component of the momentum relevant for rescaling/reversal is the
projection of $p$ onto the direction
$\hat{\delta}_{ab}=\delta_{ab}/|\delta_{ab}|$ of a vector with elements

$\delta_{ab}^{j}={2\over
m_{j}}\sum_{a^{\prime}}\operatorname{Re}\left[c_{a^{\prime}}^{*}d^{j}_{a^{\prime}a}c_{a}-c_{a^{\prime}}^{*}d^{j}_{a^{\prime}b}c_{b}\right].$
(57)

The rescaling/reversal can then be regarded as the result of a classical
particle incident on a step barrier, as in the two-level case.Mannouch and
Richardson (2023)

A key numerical advantage of MASH over conventional FSSH is that the MASH
trajectories are deterministic. This means that we can home in on each hopping
event to improve the accuracy of the time integration in a way that does not
work with stochastic hops. In the event of a hop, energy is not necessarily
conserved even after momentum rescaling, because the nuclei will have evolved
part of the time step on the ‘wrong’ surface. To resolve this, the original
step (which predicted a hop from $a$ to $b$) can be divided into two shorter
steps, one evolving on state $a$ and the other on state $b$, with the switch
time chosen as the zero of the linear approximation

$P_{b}(t)-P_{a}(t)\approx P_{b}(0)-P_{a}(0)+t\frac{P_{b}(\Delta
t)-P_{a}(\Delta t)}{\Delta t}.$ (58)

This algorithm was found to give good energy conservation for the majority of
our trajectories, even with a relatively large time step (we used a basic time
step of $\Delta t=20$ a.u. in all of the present calculations). However,
multiple hops can occasionally occur during the same time step, invalidating
the assumptions of the algorithm. We resolved this by monitoring the energy
conservation after each time step and recursively bisecting the time step
until energy was conserved to within 0.1%. If this tolerance was not reached
within 10 bisections, the trajectory was discarded (this happened for roughly
1% of the trajectories in the present calculations). Since 10 bisections
reduce the time step by a factor of 1024, this is expensive locally, but since
multiple hops in quick succession are rare events, the overall algorithm is
significantly less expensive than using such a small time step for the entire
trajectory.

## References

  * Tully (1990) J. C. Tully, “Molecular dynamics with electronic transitions,” J. Chem. Phys. 93, 1061–1071 (1990).
  * Crespo-Otero and Barbatti (2018) R. Crespo-Otero and M. Barbatti, “Recent advances and perspectives on nonadiabatic mixed quantum–classical dynamics,” Chem. Rev. 118, 7026–7068 (2018).
  * Menger  _et al._ (2018) M. F. S. J. Menger, F. Plasser, B. Mennucci, and L. González, “Surface hopping within an exciton picture. An electrostatic embedding scheme,” J. Chem. Theory Comput. 14, 6139–6148 (2018).
  * Peng  _et al._ (2022) W.-T. Peng, D. Brey, S. Giannini, D. Dell’Angelo, I. Burghardt, and J. Blumberger, “Exciton dissociation in a model organic interface: Excitonic state-based surface hopping versus multiconfigurational time-dependent Hartree,” J. Phys. Chem. Lett. 13, 7105–7112 (2022).
  * Toldo  _et al._ (2023) J. M. Toldo, M. T. do Casal, E. Ventura, S. A. do Monte, and M. Barbatti, “Surface hopping modeling of charge and energy transfer in complex environments,” Phys. Chem. Chem. Phys. 25, 8293–8316 (2023).
  * Subotnik, Ouyang, and Landry (2013) J. E. Subotnik, W. Ouyang, and B. R. Landry, “Can we derive Tully’s surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence,” J. Chem. Phys. 139, 214107 (2013).
  * Wang  _et al._ (2020) L. Wang, J. Qiu, X. Bai, and J. Xu, “Surface hopping methods for nonadiabatic dynamics in extended systems,” WIREs Comput. Mol. Sci. 10, e1435 (2020).
  * Parandekar and Tully (2006) P. V. Parandekar and J. C. Tully, “Detailed balance in Ehrenfest mixed quantum-classical dynamics,” J. Chem. Theory Comput. 2, 229–235 (2006).
  * Meyer and Miller (1979) H.-D. Meyer and W. H. Miller, “A classical analog for electronic degrees of freedom in nonadiabatic collision processes,” J. Chem. Phys. 70, 3214–3223 (1979).
  * Stock and Thoss (1997) G. Stock and M. Thoss, “Semiclassical description of nonadiabatic quantum dynamics,” Phys. Rev. Lett. 78, 578–581 (1997).
  * Runeson and Richardson (2020) J. E. Runeson and J. O. Richardson, “Generalized spin mapping for quantum-classical dynamics,” J. Chem. Phys. 152, 084110 (2020).
  * Cotton and Miller (2013) S. J. Cotton and W. H. Miller, “Symmetrical windowing for quantum states in quasi-classical trajectory simulations,” J. Phys. Chem. A 117, 7190–7194 (2013).
  * Miller and Cotton (2016) W. H. Miller and S. J. Cotton, “Classical molecular dynamics simulation of electronically non-adiabatic processes,” Faraday Discuss. 195, 9–30 (2016).
  * Mannouch and Richardson (2023) J. R. Mannouch and J. O. Richardson, “A mapping approach to surface hopping,” J. Chem. Phys. 158 (2023), 104111.
  * Kapral (2015) R. Kapral, “Quantum dynamics in open quantum-classical systems,” J. Phys. Condens. Matter 27, 073201 (2015).
  * Martens (2016) C. C. Martens, “Surface hopping by consensus,” J. Phys. Chem. Lett. 7, 2610–2615 (2016).
  * Subotnik  _et al._ (2016) J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, “Understanding the surface hopping view of electronic transitions and decoherence,” Annu. Rev. Phys. Chem. 67, 387–417 (2016).
  * Carof, Giannini, and Blumberger (2017) A. Carof, S. Giannini, and J. Blumberger, “Detailed balance, internal consistency, and energy conservation in fragment orbital-based surface hopping,” J. Chem. Phys. 147 (2017), 214113.
  * Parandekar and Tully (2005) P. V. Parandekar and J. C. Tully, “Mixed quantum-classical equilibrium,” J. Chem. Phys. 122, 094102 (2005).
  * Schmidt, Parandekar, and Tully (2008) J. R. Schmidt, P. V. Parandekar, and J. C. Tully, “Mixed quantum-classical equilibrium: Surface hopping,” J. Chem. Phys. 129, 044104 (2008).
  * Stock and Thoss (2005) G. Stock and M. Thoss, “Classical description of nonadiabatic quantum dynamics,” Adv. Chem. Phys. 131, 243–376 (2005).
  * Runeson and Richardson (2019) J. E. Runeson and J. O. Richardson, “Spin-mapping approach for nonadiabatic molecular dynamics,” J. Chem. Phys. 151, 044119 (2019).
  * Runeson  _et al._ (2022a) J. E. Runeson, J. R. Mannouch, G. Amati, M. R. Fiechter, and J. O. Richardson, “Spin-mapping methods for simulating ultrafast nonadiabatic dynamics,” Chimia 76, 582 (2022a).
  * Amati, Runeson, and Richardson (2023) G. Amati, J. E. Runeson, and J. O. Richardson, “On detailed balance in nonadiabatic dynamics: From spin spheres to equilibrium ellipsoids,” J. Chem. Phys. 158 (2023), 064113.
  * Runeson and Richardson (2021) J. E. Runeson and J. O. Richardson, “Quantum entanglement from classical trajectories,” Phys. Rev. Lett. 127, 250403 (2021).
  * Note (1) In their paper, Mannouch and Richardson used spin vectors rather than $c$ coordinates, but these contain the same information and are simply a different a notation for the same quantity.
  * Bonella and Coker (2003) S. Bonella and D. F. Coker, “Semiclassical implementation of the mapping Hamiltonian approach for nonadiabatic dynamics using focused initial distribution sampling,” J. Chem. Phys. 118, 4370–4385 (2003).
  * Müller and Stock (1999) U. Müller and G. Stock, “Flow of zero-point energy and exploration of phase space in classical simulations of quantum relaxation dynamics. II. Application to nonadiabatic processes,” J. Chem. Phys. 111, 77–88 (1999).
  * Ishizaki and Fleming (2009) A. Ishizaki and G. R. Fleming, “Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature,” Proc. Nat. Acad. Sci. 106, 17255–17260 (2009).
  * Sindhu and Jain (2022) A. Sindhu and A. Jain, “Coherence and efficient energy transfer in molecular wires: Insights from surface hopping simulations,” ChemPhysChem 23, e202200392 (2022).
  * Mannouch and Richardson (2020) J. R. Mannouch and J. O. Richardson, “A partially linearized spin-mapping approach for nonadiabatic dynamics. II. Analysis and comparison with related approaches,” J. Chem. Phys. 153, 194110 (2020).
  * Tao and Miller (2010) G. Tao and W. H. Miller, “Semiclassical description of electronic excitation population transfer in a model photosynthetic system,” J. Phys. Chem. Lett. 1, 891–894 (2010).
  * Kelly and Rhee (2011) A. Kelly and Y. M. Rhee, “Mixed quantum-classical description of excitation energy transfer in a model Fenna–Matthews–Olsen complex,” J. Phys. Chem. Lett. 2, 808–812 (2011).
  * Kim and Rhee (2014) H. W. Kim and Y. M. Rhee, “Improving long time behavior of Poisson bracket mapping equation: A non-Hamiltonian approach,” J. Chem. Phys. 140, 184106 (2014).
  * Kim, Lee, and Rhee (2014) H. W. Kim, W.-G. Lee, and Y. M. Rhee, “Improving long time behavior of Poisson bracket mapping equation: A mapping variable scaling approach,” J. Chem. Phys. 141, 124107 (2014).
  * Saller, Kelly, and Richardson (2020) M. A. C. Saller, A. Kelly, and J. O. Richardson, “Improved population operators for multi-state nonadiabatic dynamics with the mixed quantum-classical mapping approach,” Faraday Discuss. 221, 150–167 (2020).
  * Lee, Huo, and Coker (2016) M. K. Lee, P. Huo, and D. F. Coker, “Semiclassical path integral dynamics: Photosynthetic energy transfer with realistic environment interactions,” Annu. Rev. Phys. Chem. 67, 639–668 (2016).
  * Cotton and Miller (2016) S. J. Cotton and W. H. Miller, “The symmetrical quasi-classical model for electronically non-adiabatic processes applied to energy transfer dynamics in site-exciton models of light-harvesting complexes,” J. Chem. Theory Comput. 12, 983–991 (2016).
  * Berkelbach (2020) T. C. Berkelbach, “Pyrho: A Python package for reduced density matrix techniques,” https://github.com/berkelbach-group/pyrho (2020).
  * Sindhu and Jain (2023) A. Sindhu and A. Jain, “An efficient decoherence scheme for fewest switches surface hopping method,” J. Chem. Phys. 158 (2023), 10.1063/5.0143552, 154109.
  * Jain and Subotnik (2015) A. Jain and J. E. Subotnik, “Surface hopping, transition state theory, and decoherence. II. Thermal rate constants and detailed balance,” J. Chem. Phys. 143, 134107 (2015).
  * Liu and Hanna (2018) J. Liu and G. Hanna, “Efficient and deterministic propagation of mixed quantum-classical Liouville dynamics,” J. Phys. Chem. Lett. 9, 3928–3933 (2018).
  * Cao  _et al._ (2020) J. Cao, R. J. Cogdell, D. F. Coker, H.-G. Duan, J. Hauer, U. Kleinekathöfer, T. L. C. Jansen, T. Mančal, R. J. D. Miller, J. P. Ogilvie, V. I. Prokhorenko, T. Renger, H.-S. Tan, R. Tempelaar, M. Thorwart, E. Thyrhaug, S. Westenhoff, and D. Zigmantas, “Quantum biology revisited,” Sci. Adv. 6 (2020).
  * Kundu, Dani, and Makri (2022) S. Kundu, R. Dani, and N. Makri, “Tight inner ring architecture and quantum motion of nuclei enable efficient energy transfer in bacterial light harvesting,” Sci. Adv. 8, eadd0023 (2022).
  * Runeson  _et al._ (2022b) J. E. Runeson, J. E. Lawrence, J. R. Mannouch, and J. O. Richardson, “Explaining the efficiency of photosynthesis: Quantum uncertainty or classical vibrations?” J. Phys. Chem. Lett. 13, 3392–3399 (2022b).
  * Sarovar  _et al._ (2010) M. Sarovar, A. Ishizaki, G. R. Fleming, and K. B. Whaley, “Quantum entanglement in photosynthetic light-harvesting complexes,” Nat. Phys. 6, 462 (2010).
  * Lawrence  _et al._ (2019) J. E. Lawrence, T. Fletcher, L. P. Lindoy, and D. E. Manolopoulos, “On the calculation of quantum mechanical electron transfer rates,” J. Chem. Phys. 151, 114119 (2019).
  * Shushkov, Li, and Tully (2012) P. Shushkov, R. Li, and J. C. Tully, “Ring polymer molecular dynamics with surface hopping,” J. Chem. Phys. 137 (2012), 22A549.
  * Adolphs and Renger (2006) J. Adolphs and T. Renger, “How proteins trigger excitation energy transfer in the FMO complex of green sulfur bacteria,” Biophys. J. 91, 2778 – 2797 (2006).
  * Hele (2011) T. J. H. Hele,  _An Electronically Non-Adiabatic Generalization of Ring Polymer Molecular Dynamics_ , Master’s thesis, University of Oxford (2011).