# Mean-Field Ring Polymer Rates Using a Population Dividing Surface

Nathan London  [  Siyu Bu  [  Britta Johnson  [  Nandini Ananth
<EMAIL_ADDRESS>[

###### Abstract

Mean-field Ring Polymer Molecular Dynamics (MF-RPMD) offers a computationally
efficient method for the simulation of reaction rates in multi-level systems.
Previous work has established that, to model a nonadiabatic state-to-state
reaction accurately, the dividing surface must be chosen to explicitly sample
kinked ring polymer configurations where at least one bead is in a different
electronic state than the others. Building on this, we introduce a population
difference coordinate and a kink-constrained dividing surface, and we test the
accuracy of the resulting mean-field rate theory on a series of linear
vibronic coupling model systems as well as spin-boson models. We demonstrate
that this new MF-RPMD rate approach is efficient to implement and
quantitatively accurate for models over a wide range of driving forces,
coupling strengths, and temperatures.

###### keywords:

American Chemical Society, LaTeX

Cornell University]Department of Chemistry, Cornell University, Ithaca, NY
Current address: Division of Energy, Matter and Systems, School of Science and
Engineering, University of Missouri-Kansas City, Kansas City, MO Cornell
University]Department of Chemistry, Cornell University, Ithaca, NY Cornell
University]Department of Chemistry, Cornell University, Ithaca, NY Current
address: Physical Science Division, Pacific Northwest National Laboratory
Cornell University]Department of Chemistry, Cornell University, Ithaca, NY

##  1 Introduction

Nonadiabatic condensed phase reactions play a critical role in understanding
reaction mechanisms for a range interesting systems including proton coupled
electron transfer in biological systems and charge transfer in energetic
materials. 1, 2, 3, 4, 5, 6 Despite the development of numerous nonadiabatic
dynamic methods for the characterization of reaction rates and mechanisms,
atomistic simulations of multi-level systems that incorporate nuclear quantum
effects and vibronic couplings remain largely out of reach.

Perhaps the most ubiquitous and successful nonadiabatic rate theory is the
Marcus theory for electron transfer that assumes a classical solvent. 7 Going
beyond Marcus theory, Fermi’s Golden Rule (FGR) rate theory incorporates
nuclear quantum effects, and can be computed exactly for parabolic potentials.
In certain regimes, FGR rates can be approximated using Wolynes rate theory, 8
and limitations of the Wolynes approach have been addressed by recent work,
extending the methods applicability to the Marcus inverted regime and
improving its accuracy for high temperatures and anharmonic potentials. 9, 10,
11, 12, 13 Despite these advances, these theories are confined to the weak
perturbation limit where FGR is applicable, although an interpolation strategy
has been proposed to move between the golden-rule limit theories and adiabatic
rate theory approaches. 14

Dynamical methods developed to study the mechanisms of nonadiabatic processes
and to compute rate constants include exact methods like MultiConfiguration
Time-Dependent Hartree (MCTDH) 15, 16 real-time path-integral based
approaches, 17, 18, 19 and Heirarchical Equations of Motion, 20, 21, 22 in
addition to approximate methods like Redfield theory, 23, 24 Fewest Switches
Surface Hopping and other mixed quantum-classical dynamics methods, 25, 26,
27, 28, 29, 30, 31 and Semiclassical theories. 32, 33, 34, 35, 36, 37, 38, 39
Unfortunately, the majority of these approaches are not suitable for large
scale system simulations, with a notable exception being systems that can be
modeled by spin-boson Hamiltonians where transitions between electronic states
are mediated by coupling to a nuclear bath. Approximate imaginary-time path
integral methods like Centroid Molecular Dynamics (CMD) 40, 41 and Ring
Polymer Molecular Dynamics (RPMD) 42 show great promise in moving beyond the
spin-boson models to the more complex bilinear coupling Hamiltonian models and
even atomistic simulations. 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
55, 56, 57, 58 . Both CMD and RPMD conserve the quantum Boltzmann distribution
and capture nuclear quantum effects like tunneling using an ensemble of
classical trajectories in an extended phase space, making them efficient for
the simulation of systems with many degrees of freedom. Nonadiabatic
extensions of CMD have been explored in the context of the spin-boson model,
53 while RPMD has been used for atomistic simulations of nonadiabatic electron
transfer and proton-coupled electron transfer 50, 52 However, it has been
shown that RPMD cannot capture state-to-state dynamics when working with a
two-level or, in general, $K$-level system Hamiltonian. 50 This limitation has
motivated the development of several multi-state imaginary-time path integral
based methods including mean-field (MF)-RPMD 59, 60, 61, nonadiabatic RPMD 62,
kinetically constrained (KC)-RPMD 63, 64, 65, coherent-state RPMD 66, and
mapping-variable (MV)-RPMD. 67, 68, 69, 70 Other approaches have also been
developed that combine RPMD with adhoc nonadiabatic dynamics methods like
surface hopping. 71, 72, 73, 74

Of these methods, only MV-RPMD and MF-RPMD retain the key RPMD feature of
conserving the quantum Boltzmann distribution, and MF-RPMD is unique in that
it does not introduce any explicit electronic variables, instead describing
nuclear dynamics on an effective state-averaged potential energy surface. This
makes MF-RPMD more approximate but significantly more efficient than MV-RPMD.
Initial efforts to compute MF-RPMD rates used the nuclear centroid reaction
coordinate and defined the dividing surface as configurations where the two
electronic states were degenerate. 60 It was shown that the resulting rate
expression significantly overestimated nonadiabatic reaction rates. 60
Subsequent work has established that this error can be attributed to the use
of a dividing surface that does not require the formation of kinked ring
polymer configurations (where at least one bead is in a different electronic
state than the others), an essential intermediate step in state-to-state
reactions. 64, 61 Further, it has been established that this error can be
mitigated by introducing a doubly-constrained dividing surface that requires
formation of kinked configurations, however the resulting MF-RPMD rate
expression proved challenging to implement. 61

Here we introduce an improved reaction coordinate that tracks the difference
in electronic state populations, and a modified reactive flux definition that
accounts for the formation of kinked ring polymer configurations. We
demonstrate that this new coordinate and reactive flux definition enable
accurate calculations of nonadiabatic and adiabatic reaction rates for a
series of bilinear coupling model systems, captures nonadiabatic electron
transfer rates across a wide range of driving forces, and successfully
incorporates nuclear quantum effects by reproducing FGR rates over a range of
temperatures in models where the driving force places them in the Marcus
normal regime. In the Marcus inverted regime, we show MF-RPMD reproduces
Marcus theory rates rather than FGR rates. This failure to account for nuclear
tunneling can be traced to the use of a Boltzmann weighted reaction coordinate
that effectively limits the the ring polymer to sample only nuclear
configurations where the two electronic states are degenerate.

##  2 Methods

###  2.1 Path Integrals for Multi-level Systems

For a general $K$-level system with $d$ nuclear degrees of freedom, the
diabatic Hamiltonian is

$\hat{H}=\sum_{j=1}^{d}\frac{\hat{P}_{j}^{2}}{2M_{j}}+\sum_{n,m=1}^{K}|n\rangle
V_{nm}(\hat{R})\langle m|,$ (1)

where $\hat{R}$ and $\hat{P}$ are $d$-dimensional nuclear position and
momentum vector operators, respectively, $M_{j}$ is the nuclear mass,
$\\{|n\rangle\\}$ is the set of diabatic states, and $V_{nm}(\hat{R})$ is the
diabatic potential energy operator. The quantum canonical partition function
is obtained by evaluating the trace in the basis of nuclear positions and
electronic states, and discretized via repeated insertion of $N$ copies of
identity in the form of completeness relations to obtain

$\displaystyle Z$ $\displaystyle=$
$\displaystyle\mathrm{Tr}[e^{-\beta\hat{H}}]$ (2) $\displaystyle=$
$\displaystyle\int
d\\{R_{\alpha}\\}\sum_{\\{n_{\alpha}\\}=1}^{K}\prod_{\alpha=1}^{N}\langle
R_{\alpha},n_{\alpha}|e^{-\beta_{N}\hat{H}}|R_{\alpha+1},n_{\alpha+1}\rangle,$

where $\beta_{N}=1/(Nk_{B}T)$, T is temperature, N is the number of imaginary
time slices (or beads), and $R_{\alpha}$, $n_{\alpha}$ refer to the nuclear
position and electronic state, respectively, of the $\alpha^{\text{th}}$ bead.
In Eq. 2, the shorthand, $\int d\left\\{R_{\alpha}\right\\}=\int dR_{1}\int
dR_{2}..\int dR_{N}$ and
$\sum_{\\{n_{\alpha}\\}=1}=\sum_{n_{1}}\sum_{n_{2}}\ldots\sum_{n_{N}}$, is
used to represent the multi-dimensional integral over nuclear coordinates and
summation over electronic states. The matrix elements in Eq. 2 are evaluated
using the Trotter (high temperature) approximation to obtain 75, 76,

$\displaystyle
Z\propto\lim_{N\rightarrow\infty}\int\\{dR_{\alpha}\\}e^{-\beta_{N}V_{N}(\\{R_{\alpha}\\})}\text{Tr}[\Gamma],$
(3)

where

$V_{N}=\sum_{j=1}^{d}\sum_{\alpha=1}^{N}\left[\frac{M_{j}}{2\beta_{N}^{2}}(R_{j,\alpha}-R_{j,\alpha+1})^{2}\right],$
(4) $\Gamma=\prod_{\alpha=1}^{N}\mathbf{M}(R_{\alpha}),$ (5)

and $\mathbf{M}$ is the $K\times K$-dimensional matrix

$M_{nm}(R_{\alpha})=\begin{cases}e^{-\beta_{N}V_{nn}(R_{\alpha})}&n=m\\\
-\beta_{N}V_{nm}(R_{\alpha})e^{-\beta_{N}V_{nn}(R_{\alpha})}.&n\neq
m\end{cases}$ (6)

The trace in Eq. 3 can be moved into the exponential, and $N$ normalized
nuclear momentum gaussians can be introduced to express the canonical
partition function as a phase space integral,

$Z\propto\lim_{N\rightarrow\infty}\int\\{dR_{\alpha}\\}\int\\{dP_{\alpha}\\}e^{-\beta_{N}H_{N}(\\{R_{\alpha}\\},\\{P_{\alpha}\\})},$
(7)

where the Mean-Field Ring Polymer Hamiltonian is

$H_{N}=\sum_{j=1}^{d}\sum_{\alpha=1}^{N}\left[\frac{M_{j}}{2\beta_{N}^{2}}(R_{j,\alpha}-R_{j,\alpha+1})^{2}+\frac{P_{j,\alpha}^{2}}{2M_{j}}\right]-\frac{1}{\beta_{N}}\ln{(\text{Tr}[\Gamma])}.$
(8)

In Eq. 8, the Hamiltonian is a function of nuclear phase space variables only,
with an effective mean-field potential,
$V_{\text{MF}}=-\frac{1}{\beta_{N}}\ln{(\text{Tr}[\Gamma])}$, obtained by
tracing over all possible electronic state configurations for the ring
polymer. During a state-to-state reaction, the contributions to the effective
potential from ring polymer configurations where all the beads are in a single
electronic state and configurations where the beads are in different
electronic states will vary as the system moves from the reactant to product
state as shown in Fig. 1.

Figure 1: Cartoon showing the progress of a multi-level system reaction with
state-space path integrals. In the reactant and product states, ring polymer
configurations where all beads are in a single electronic state are the
dominant contributions to the effective mean-field potental. At the transition
state, when the system is moving from one state to the other, we expect to see
a significant contribution from multi-electronic state or kinked ring polymer
configurations. 77

###  2.2 The MF-RPMD Approximation to Reaction Rates

The MF Hamiltonian in Eq. 8 can be used to generate molecular dynamics
trajectories that sample equilibrium configurations of the $K$-level system at
finite temperature as is done in standard Path Integral Molecular Dynamic
simulations. MF-RPMD approximates the Kubo-transformed real-time quantum
correlation function by sampling trajectory initial conditions from an exact
quantum canonical ensemble and time-evolving them under the mean-field
Hamiltonian. By construction, MF-RPMD trajectories preserve the quantum
Boltzmann distribution.

The MF-RPMD approximation to the flux-side correlation function, in keeping
with the RPMD approximation, 78 can be written as,

$C^{\text{MF}}_{\text{fs}}(t)=\frac{\langle\delta(\xi_{0}-\xi^{\ddagger})\dot{\xi_{0}}h(\xi_{t}-\xi^{\ddagger})\rangle}{\langle
h(\xi^{\ddagger}-\xi_{0})\rangle},$ (9)

where $\langle\dots\rangle$ indicates a canonical ensemble average. The
reaction coordinate, $\xi$, in Eq. 9 is a function of the ring polymer nuclear
configurations $\\{R\\}$, with $\xi_{0}$ being its value at time $t=0$,
$\xi_{t}$ the value at any later time $t$, and $\xi^{\ddagger}$ defines the
dividing surface. For reactions involving transitions between electronic
states, although there are no explicit electronic state variables in the MF-
RPMD Hamiltonian, previous work suggests that it is essential that the
definition of the reactive flux properly accounts for the formation of
‘kinked’ ring polymer configurations as shown in Fig 1. 77, 63, 61
Unfortunately, for nonadiabatic processes, sampling the dividing surface
$\xi=\xi^{\ddagger}$ in Eq. 9 with the mean-field Hamiltonian yields an
ensemble of ring polymer configurations where all beads are in a single
electronic state rather than kinked configurations.  The origin of this
difficulty is well understood: the probability of kink formation is
proportional to the square of the diabatic coupling, making this a rare event
in the nonadiabatic regime where coupling values are very small. In the
language of the mean-field effective potential,

$\displaystyle\Gamma=\Gamma_{1}+\Gamma_{2}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{\mathrm{k}}},$
(10)

where $\Gamma_{i}$, $i=1,2$, represent the potential contribution from
configurations where all ring polymer beads are in the same electronic state,

$\displaystyle\Gamma_{i}=\prod_{\alpha=1}^{N}\left[\mathbf{M}(R_{\alpha})\mathcal{P}_{i}\right],$
(11)

where $\mathcal{P}_{i}$ is the projection operator matrix corresponding to the
$i^{\text{th}}$ electronic state, and the matrix elements of $\mathbf{M}$ are
defined in Eq. 6. The contribution to the effective potential from kinked ring
polymer configurations,
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Gamma_{\mathrm{k}}}$,
has a coupling-square dependence that ensures that it is _not_ the dominant
contribution to the $\Gamma$ function at any nuclear configuration.

To sample _only_ kinked configurations at the dividing surface, it is
necessary to introduce an adhoc modification to the dividing surface that re-
weights the sampling,

$\displaystyle\delta(\xi-\xi^{\ddagger})\rightarrow\delta(\xi-\xi^{\ddagger})\frac{\text{Tr}\left[{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>(12)

Using the modified dividing surface, the MF-RPMD rate can be expressed as,

$k=\lim_{t\rightarrow t_{\text{plateau}}}\frac{\int
d\\{R_{\alpha},P_{\alpha}\\}e^{-\beta_{N}H_{N}^{0}}\text{Tr}\left[\Gamma_{\mathrm{k}}\right]\delta(\xi_{0}-\xi^{\ddagger})\dot{\xi_{0}}h(\xi_{t}-\xi^{\ddagger})}{\int
d\\{R_{\alpha},P_{\alpha}\\}e^{-\beta_{N}H_{N}^{0}}\text{Tr}\left[\Gamma_{1}\right]},$
(13)

where $H_{N}^{0}$ is the free ring polymer Hamiltonian that includes only the
harmonic spring term and the kinetic energy terms of Eq. 8 and the presence of
$\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\Gamma_{\mathrm{k}}}$
ensures that the effective potential samples only kinked configurations. The
time-evolved reaction coordinate, $\xi_{t}$, is obtained by integrating the
classical equations of motion due to the MF-RPMD Hamiltonian in Eq. 8. The
denominator in Eq. 13 is the reactant partition function, and is calculated by
sampling configurations projected onto the reactant (here, state 1).

We now introduce a new reaction coordinate, a bead-averaged population
difference estimator that was previously employed in the context of MF-RP
instanton calculations, 59, 79 and is exact at equilibrium,

$\Delta
P=\frac{1}{N}\sum_{\alpha=1}^{N}\frac{e^{-\beta_{N}V_{22}(R_{\alpha})}-e^{-\beta_{N}V_{11}(R_{\alpha})}}{e^{-\beta_{N}V_{11}(R_{\alpha})}+e^{-\beta_{N}V_{22}(R_{\alpha})}}.$
(14)

This coordinate takes negative values in the reactant region, positive values
in the product region, and is zero at the dividing surface. The corresponding
MF-RPMD rate expression is then,

$k=\lim_{t\rightarrow t_{\text{plateau}}}\frac{\int
d\\{R_{\alpha},P_{\alpha}\\}e^{-\beta_{N}H_{N}^{0}}\text{Tr}\left[\Gamma_{\mathrm{k}}\right]\delta(\Delta
P_{0})\dot{\Delta P_{0}}h(\Delta P_{t})}{\int
d\\{R_{\alpha},P_{\alpha}\\}e^{-\beta_{N}H_{N}^{0}}\text{Tr}\left[\Gamma_{1}\right]},$
(15)

The population coordinate and dividing surface introduced here are
significantly different from the coordinate previously introduced by one of
us. 61 In that work, the dividing surface was defined as nuclear
configurations where half the ring polymer beads were in one state and the
other half were in a different state. In addition, the reaction coordinate
measured population differences but with the so-called Boltzmann estimator, 70
is given by

$\displaystyle\chi=\frac{\text{Tr}\left[\prod_{\alpha=1}^{N}{\bf
M}(R_{\alpha})\left|2\right\rangle\left\langle
2\right|\right]-\text{Tr}\left[\prod_{\alpha=1}^{N}{\bf
M}(R_{\alpha})\left|1\right\rangle\left\langle
1\right|\right]}{\text{Tr}\left[\prod_{\alpha=1}^{N}{\bf
M}(R_{\alpha})\left|2\right\rangle\left\langle
2\right|\right]+\text{Tr}\left[\prod_{\alpha=1}^{N}{\bf
M}(R_{\alpha})\left|1\right\rangle\left\langle 1\right|\right]}$ (16)

While both population estimators are exact at equilibrium in the mean-field
path integral framework, the coordinate used here, $\Delta P$ as defined Eq.
14, has been previously used to track changes in electronic state populations
along the mean-field ring polymer instantons. 79, 59 There are two significant
advantages to the present $\Delta P$ coordinate. First, this coordinate is a
bead-average population difference coordinate ensuring that the time zero
correlation function correctly corresponds to a Kubo-transform real-time
thermal correlation function. Second, $\Delta P$ can be easily computed in
large scale simulations where the donor and acceptor state potentials are
known ($V_{11}$ and $V_{22}$), and does not require the matrix multiplication
steps necessary to construct the $\chi$ coordinate. As a final point of
comparison, we note that the velocity corresponding to $\chi$ was computed
numerically using finite difference, 61 whereas the analytic time derivative
of $\Delta P$ is easily computed.

Our new rate expression defined in Eq. (15) has some parallels with the rate
formalism of KC-RPMD. 64 Both methods emphasize the formation of kinked RP
configurations to accurately described nonadiabatic reaction rates, but differ
in the details. In KC-RPMD, the reaction coordinate is defined by an auxiliary
dynamical variable that reports on the formation of kink-pairs; the
introduction of this variable results in a modified Boltzmann distribution
that the KC-RPMD Hamiltonian conserves. The MF-RPMD rate presented here does
not introduce any additional dynamical variables, with our reaction coordinate
depending only on the RP bead positions and the diabatic potential energy
surfaces. Further, in KC-RPMD, the restriction that nonadiabatic transitions
should occur at configurations where the diabatic states are near-degenerate
is enforced by the use of a penalty function on kinked configurations in the
effective potential in the form of a Gaussian of the scaled energy difference
of the diabatic surfaces. In the present work, we ensure kink formation near
the diabatic state crossing through our modified dividing surface.

###  2.3 Reference Rate Theories

We compare nonadiabatic MF-RPMD rates computed here against Marcus theory and
FGR rates. The Marcus theory rate for a nonadiabatic process with a classical
solvent is 1,

$k_{\text{MT}}=\frac{2\pi}{\hbar}{|V_{nm}|}^{2}\sqrt{\frac{\beta}{4\pi\lambda}}e^{-\beta\left(\lambda-\varepsilon\right)^{2}/4\lambda},$
(17)

where $\lambda$ is the reorganization energy, $\varepsilon$ is the driving
force, and $V_{nm}$ is the diabatic coupling between the reactant and product
electronic states.

FGR rate theory for a nonadiabatic reaction using the flux correlation
formalism is 80,

$k_{\mathrm{FGR}}=\frac{1}{Z_{0}}\frac{\Delta^{2}}{\hbar^{2}}\int_{-\infty}^{\infty}dt\
C^{\tau}(t),$ (18)

where

$C^{\tau}(t)=\mathrm{Tr}\left[\mathrm{e}^{i\hat{H}_{0}(t+i(\beta\hbar-\tau))/\hbar}\mathrm{e}^{-i\hat{H}_{1}(t-i\tau)/\hbar}\right],$
(19)

$Z_{0}=\mathrm{Tr}\left[\mathrm{e}^{-\beta\hat{H}_{0}}\right]$ is the reactant
partition function. $\hat{H}_{0}$ and $\hat{H}_{1}$ are the Hamiltonians
representing the reactant and product states respectively. $\tau$ can be any
real number as the rate is independent of its choice 80 but is, typically,
chosen to be in the range $[0,\beta\hbar]$. For the bilinear coupling model
Hamiltonians employed here, the FGR rate can be obtained analytically; 81
Details of this calculation are provided in the Supporting Information.

In the adiabatic limit, the electron transfer rate constant for a system with
quantum solvent can be expressed using Kramers theory, 82

$k_{\mathrm{KT}}=\left(\sqrt{1+\left(\frac{\gamma}{2\omega_{\mathrm{b}}}\right)^{2}}-\frac{\gamma}{2\omega_{\mathrm{b}}}\right)\frac{\omega_{\mathrm{s}}}{2\pi}\mathrm{e}^{-\beta
G_{\mathrm{cl}}^{\ddagger}},$ (20)

where $\omega_{\mathrm{s}}$ is the frequency of the harmonic potential,
$\omega_{\mathrm{b}}$ is the frequency at the barrier top,
$\gamma=\eta/M_{\mathrm{s}}$, $\eta$ is the strength of the coupling to a
dissipative bath, $M_{\mathrm{s}}$ is the solvent mass, and
$G_{\mathrm{cl}}^{\ddagger}$ is the free energy barrier computed for a
classical solvent.

###  2.4 Model Systems

We calculate rates for a series of bilinear vibronic coupling models of
condensed phase electron transfer with potential

$V(\hat{R})=V_{\mathrm{S}}(\hat{s})+V_{\mathrm{B}}(\hat{R})$ (21)

where the configuration vector $\hat{R}=\\{\hat{s},\hat{Q}\\}$ represents the
solvent polarization coordinate, $s$, and the bath coordinates, $Q$. The
diabatic potential energy matrix is

$V_{\mathrm{S}}(\hat{s})=\left(\begin{array}[]{cc}A\hat{s}^{2}+B\hat{s}+\varepsilon&\Delta\\\
\Delta&A\hat{s}^{2}-B\hat{s}\end{array}\right),$ (22)

where $\varepsilon$ is the driving force of the reaction and $\Delta$ is the
constant diabatic coupling. The solvent coordinate is linearly coupled to a
thermal bath of $f$ harmonic oscillators,

$V_{\mathrm{B}}(\hat{R})=\sum_{j=1}^{f}\left[\frac{1}{2}M_{\mathrm{B}}\omega_{j}^{2}\left(\hat{Q}_{j}-\frac{c_{j}\hat{s}}{M_{\mathrm{B}}\omega_{j}^{2}}\right)^{2}\right]$
(23)

where $M_{B}$ is the bath mass. The bath is described by an Ohmic spectral
density

$J(\omega)=\eta\omega e^{-\frac{\omega}{\omega_{c}}}$ (24)

where $\omega_{c}$ is the cutoff frequency and $\eta$ is the dimensionless
friction coefficient. The spectral density is discretized into $f$ oscillators
43

$\omega_{j}=-\omega_{\mathrm{c}}\ln\left(\frac{j-0.5}{f}\right),$ (25)

with coupling strengths

$c_{j}=\omega_{j}\left(\frac{2\eta
M_{\mathrm{B}}\omega_{\mathrm{c}}}{f\pi}\right)^{1/2}.$ (26)

The parameters for this model are shown in Table 1.

Table 1: Linear vibronic coupling model parameters in atomic units unless otherwise indicated. Parameters | Value
---|---
$A$ | $4.772\times 10^{-3}$
$B$ | $2.288\times 10^{-2}$
$\varepsilon$ | $0.0-0.2366$
$\Delta$ |  $6.69\times 10^{-7}$ \- $1.20\times 10^{-2}$
$M_{\mathrm{S}}$ | $1836.0$
$M_{\mathrm{B}}$ | $1836.0$
$f$ | $12$
$\omega_{\mathrm{c}}$ | $2.28\times 10^{-3}$
$\eta/M_{B}\omega_{c}$ | 1.0
$T$ |  $150$ \- $300$ K

We further investigate a set of spin-boson models previously benchmarked using
Heirarchy Equations of Motion (HEOM) against RPMD and Wolynes theory. 74 These
models allow us to investigate the performance of MF-RPMD in the overdamped
and underdamped regimes for both high frequency (B2 and B4) and low frequency
(B1 and B3) systems. Parameters for all four models are shown in Table 2, with
the shared parameters reported in Table 3.

Table 2: Parameters (in atomic units) used for the diabatic potential in Eq. (21) that defines the spin-boson models B1-B4. Parameters | B1 | B2 | B3 | B4
---|---|---|---|---
$A$ | 0.125 | 8 | 0.125 | 8
$B$ | 2.739 | 21.909 | 2.739 | 21.909
$\eta/M_{B}\omega_{c}$ | 32 | 32 | 1 | 1
$\omega_{\mathrm{c}}$ | 0.5 | 4 | 0.5 | 4
Table 3: Shared parameters for the spin-boson models B1-B4 in atomic units. Parameters | Value
---|---
$\varepsilon$ | $15$
$\Delta$ | 0.1585, 1, 2.5119
$M_{\mathrm{S}}$ | 1.0
$M_{\mathrm{B}}$ | 1.0
$f$ | $12$
$\beta$ | 1

###  2.5 Simulation Details

The numerator of the rate in Eq. 13 requires constrained sampling to the
dividing surface while the denominator requires sampling the reactant region.
We ensure sampling of all relevant configuration space using a windowing
approach that introduces identity in the form of an integral over all solvent
ring polymer centroid configurations. The resulting rate expression is

$k=\lim_{t\rightarrow t_{\text{plateau}}}\frac{\int
ds^{\prime}\langle\text{Tr}[\Gamma_{\mathrm{k}}]\delta(\Delta
P_{0})\dot{\Delta P}_{0}h(\Delta P_{t})\rangle_{w}}{\int
ds^{\prime}\langle\text{Tr}[\Gamma_{1}]\rangle_{w}}$ (27)

where $\langle...\rangle_{w}$ is used to indicate a phase space ensemble
average over the nuclear bead configurations obtained by importance sampling
from the constrained distribution

$w=e^{-\beta_{N}H_{N}^{0}\left(\left\\{R_{\alpha},P_{\alpha}\right\\}\right)-k_{\mathrm{s}}(\bar{s}-s^{\prime})^{2}},$
(28)

where $\bar{s}$ is the solvent ring polymer centroid. The integral in Eq. 27
is evaluated performed by ensemble averaging for a series of windows centered
at different values of $s^{\prime}$, and then using the trapezoid rule. A
standard Metropolis algorithm is used to sample free ring polymer
configurations from
$e^{-\beta_{N}H_{N}^{0}\left(\left\\{R_{\alpha},P_{\alpha}\right\\}\right)}$
subject to a harmonic constraint with force constant $k=1000$ a.u. on the
centroid. We use 150 evenly spaced windows over a range of $s^{\prime}$
values. The initial value of $s^{\prime}$ for each model is given in Table S1
in the Supporting Information.

All MF-RPMD simulations presented here are converged with respect to bead
number, as reported in the results section, convergence numbers range from
$N=16-32$. Initial conditions are sampled from $5\times 10^{4}$ decorrelated
Monte Carlo steps for each window, and time-evolved under the mean-field
Hamiltonian in Eq. 8 using a fourth-order Adams-Bashforth-Moulton predictor-
corrector integrator with a time step of 0.05 a.u. for models I-VIII and
0.0005 a.u. for models B1-4. The flux-side correlation function plateau in
time is achieved between 500 a.u. and 3000 a.u. for models I-VIII and between
5 a.u. to 300 a.u. for the B1-4 models.

The delta function at time zero in Eq. 27 is approximated with the Gaussian,

$\delta(\Delta P)\approx\frac{1}{\sigma\sqrt{2\pi}}\
\mathrm{e}^{-\frac{1}{2}\left(\frac{\Delta P}{\sigma}\right)^{2}},$ (29)

where $\sigma$ is treated a parameter controlling the width of the Gaussian.
Results reported here employ $\sigma=0.01-0.001$; values that are sufficiently
small to mimic a delta function, but sufficiently large to ensure meaningful
statistics in terms of number of trajectories and number of windows
contributing to the overall rate integral. We ensure that the rates reported
are consistent across at least two smaller $\sigma$ values for each model. The
analytical gradient of the population difference coordinate is used to obtain
the initial velocity along the reaction coordinate in Eq. 9.

The flux-side correlation function in Eq. 13 assumes that the velocity is
‘positive’ when the system is moving from reactant to product. In the inverted
regime, the analytic time-derivative of the population difference coordinate
can be negative for trajectories starting on the reactant side of the minimum,
and this can yield a reaction rate that is negative. We correct this by
including a sign function that changes the sign of an individual bead’s
contribution to ensure that the bead-average velocity has a positive sign when
the trajectory is moving towards product formation,

$\displaystyle\dot{\Delta P}$
$\displaystyle=\frac{1}{N}\sum_{\alpha}\dot{\Delta P}_{\alpha}$ (30)
$\displaystyle\dot{\Delta P}_{\alpha}$ $\displaystyle=-\text{sgn}(\dot{\Delta
P}_{\alpha})$ (31)

##  3 Results and Discussion

We start with a study of MF-RPMD rates for nonadiabatic reactions with
different driving forces at temperature $T=300$ K and with the diabatic
coupling $\Delta=6.69\times 10^{-7}$. As shown in Fig. 2, we find MF-RPMD
rates in quantitative agreement with FGR rates for Models I-IV that are in the
Marcus normal regime. The results shown in the figure are also tabulated in
Table 4, with all results calculated using $\sigma=0.01$ in Eq. 29. Note that
we also calculate MF-RPMD rates using a centroid coordinate for normal regime
models I, II, and IV, but do not plot the results since they are visually
indistinguishable from the population coordinate results shown here. The rate
expression with the centroid coordinate and the resulting MF-RPMD rates are
provided in the SI.

Figure 2: Nonadiabatic rates for bilinear coupling models I-VIII. MF-RPMD rates with are shown with red circles. FGR rates (triangles and solid line) and the Marcus theory rates (squares and dashed line) are also shown. Table 4: Nonadiabatic reaction rates for models I-VIII Model | $\varepsilon$ |  log($k_{\mathrm{MF}}$) |  log($k_{\mathrm{FGR}}$) |  log($k_{\mathrm{MT}}$)
---|---|---|---|---
I | 0.00 | -21.42 | -21.37 | -22.65
II | 0.0158 | -18.09 | -18.08 | -19.30
III | 0.0296 | -15.60 | -15.73 | -16.79
IV | 0.0586 | -11.70 | -12.26 | -12.83
V | 0.1186 | -9.97 | -10.26 | -10.19
VI | 0.1776 | -14.38 | -13.22 | -14.91
VII | 0.2071 | -19.35 | -16.07 | -19.99
VIII | 0.2366 | -26.10 | -19.65 | -26.89

In the inverted regime, for models V-VIII, the MF-RPMD rate is in agreement
with the Marcus Theory rate, despite the fact that MF-RPMD quantizes the
solvent coordinate. We find that our definition of the dividing surface at
$\Delta P=0$ is responsible for this discrepancy. As the population of each
state is defined using the Boltzmann weight of that state, the state with the
lower energy will have a higher population. In the normal regime, $\Delta P=0$
can be achieved by tunneling configurations where some beads are the left of
the crossing with $\Delta P<0$ and some are to the right with $\Delta P>0$
resulting in a bead-average value of zero. In the inverted regime, this
remains true, however, tunneling configurations where all beads are to the
right of the crossing fail to meet the $\Delta P=0$ criterion as shown in the
cartoon Fig. 3. This result highlights one of the limitations of using MF-
RPMD: the lack of an explicit electronic variable means that constraining to
kinked configurations places potentially unphysical constraints on the nuclear
configurations. In contrast, although KC-RPMD also employs a kink-constrained
dividing surface, the constraint applies only to the auxiliary electronic
variable introduced in that framework, allowing the method to capture nuclear
tunneling in the Marcus inverted regime.

Figure 3: Cartoon of RP configurations in the inverted regime showing kinked
configurations where some beads are in the reactant state and others are in
the product state. Although configurations in (a) are quantum mechanically
valid and necessary to capture solvent tunneling, in our MF-RPMD simulation
they correspond to $\Delta P\neq 0$. Configurations in (b) do obey the $\Delta
P=0$ criterion. For this reason, the MF-RPMD rates reported here fail to
reproduce FGR rates in the inverted regime, instead yielding Marcus theory
rates.

One of the strengths of MF-RPMD is the ability to accurately compute both
adiabatic and nonadiabatic reaction rates within a single framework. We
demonstrate this using model I, a normal regime model at $T=300$ K, and we
plot the results calculated using $\sigma=0.01$ in Eq. 29 against the Marcus
and FGR rates in Fig. 4, and tabulate them in Table 5. We find that MF-RPMD
rates move seamlessly from good agreement with FGR rates in the nonadiabatic
cases to good agreement with Kramer’s theory rates in the adiabatic regime.

Figure 4: Reactions rates for model I with coupling strength modified to tune the system with varying coupling. MF-RPMD rates with the new population dividing surface are shown as red circles and the error bars are within the symbol size. The FGR rate (triangles and solid line) and the Kramers theory rate (squares and dashed line) are also shown Table 5: Reaction rates for model I with variable coupling $\Delta$ | $N$ |  log($k_{\mathrm{MF}}$) |  log($k_{\mathrm{FGR}}$) |  log($k_{\mathrm{KT}}$)
---|---|---|---|---
$6.69\times 10^{-7}$ | 16 | -21.42 | -21.37 | -15.57
$3.16\times 10^{-6}$ | 16 | -20.07 | -20.02 | -15.58
$3.16\times 10^{-5}$ | 16 | -18.06 | -18.02 | -15.58
$5.01\times 10^{-4}$ | 32 | -15.68 | -15.62 | -15.54
$2.00\times 10^{-3}$ | 32 | -14.42 | -14.42 | -15.02
$7.94\times 10^{-3}$ | 32 | -12.27 | -13.22 | -12.60
$1.20\times 10^{-2}$ | 32 | -11.02 | -12.86 | -11.11

Further, we examine the effect of temperature on the rate constant. For this
we choose a model in the normal regime, model II, in the weak coupling limit
$\Delta=6.69\times 10^{-7}$. For these simulations we use a value of
$\sigma=0.001$ in Eq. 29. The results are shown in Fig. 5 and tabulated in
Table 6. Consistent with our other results in the normal regime, we achieve
quantitative agreement with the rates from Fermi’s golden rule. As expected,
we do have to increase the number of beads from 16 to 24 to get converged
results at lower temperatures. At very high temperatures, as suggested by
earlier studies, 83 the present rate expression is expected to perform poorly
given the adhoc nature of the dividing surface introduced. Indeed, we find
this is the case at $T=5000K$ comparing the 8-bead MF-RPMD rate (1-bead is the
zero-coupling limit) against the corresponding FGR rate, we find the rate
differs by almost an order of magnitude, $k_{\text{FGR}}=5.36\times 10^{-12}$
a.u. and $k_{\text{MF}}=1.56\times 10^{-13}$ a.u.

Figure 5: ET rates as function of temperature for model II with $\varepsilon=0.0158$. MF-RPMD rates with the new population dividing surface are shown as red circles. The Fermi’s golden rule (triangles and solid line) and the Kramers theory (squares and dashed line) rates are also shown Table 6: MF-RPMD rates for the linear vibronic coupling model with varying temperature T (K) | $N$ |  log($k_{\mathrm{MF}}$) |  log($k_{\mathrm{FGR}}$) |  log($k_{\mathrm{MT}}$)
---|---|---|---|---
150 | 24 | -21.84 | -22.07 | -28.33
200 | 24 | -20.46 | -20.48 | -23.80
250 | 16 | -19.15 | -19.16 | -21.09
300 | 16 | -18.07 | -18.08 | -19.30

Finally, we present the MF-RPMD rates for models B1-B4. The results are shown
in Fig. 6 and tabulated in Table 7. Overall, we find relatively good agreement
when against FGR and Kramer’s rate theory results as the systems move from the
nonadiabatic to the adiabatic regime. The most significant error is for the
case of low frequency and strong system-bath coupling where we find the
nonadiabatic MF-RPMD rates significantly underestimate the FGR rates. This
appears to be in keeping with the limitations of the current rate expression
that does not correctly capture classical-limit behavior.

Figure 6: MF-RPMD reaction rates for the spin-boson models B1-B4 as a function of coupling strength. The MF-RPMD rates (red circles) are compared against with The Fermi’s golden rule (triangles and solid line) and the Kramers theory (squares and dashed line) Table 7: MF-RPMD rates for the spin-boson models B1-B4 Model |  log($\beta\Delta$) |  log($k_{\mathrm{MF}}$) |  log($k_{\mathrm{FGR}}$) |  log($k_{\mathrm{Kramer}})$
---|---|---|---|---
B1 | -0.8 | -7.03 | -5.88 | -5.59
0 | -5.58 | -4.28 | -5.77
0.4 | -4.80 | -3.48 | -5.43
B2 | -0.8 | -6.00 | -5.62 | -3.78
0 | -4.44 | -4.02 | -3.96
0.4 | -3.81 | -3.22 | -3.62
B3 | -0.8 | -6.40 | -5.87 | -5.18
0 | -4.78 | -4.27 | -5.03
0.4 | -3.89 | -3.47 | -4.50
B4 | -0.8 | -4.68 | -4.70 | -3.37
0 | -3.29 | -3.10 | -3.22
0.4 | -2.58 | -2.30 | -2.70

##  4 Conclusions

We introduce a new population difference reaction coordinate and a kink-
constrained dividing surface for the calculation of rates using MF-RPMD. This
new rate formulation represents a simplification of the previously introduced
MF-RPMD rate expression for multi-level systems with an implementation that is
designed to enable atomistic system simulations. We have shown that the MF-
RPMD rate expression introduced here is accurate over a wide range of
electronic coupling strengths, captures nonadiabatic rates accurately over a
range of driving forces, and is able to capture nuclear tunneling effects at
low temperatures.

Future directions include identifying modifications to the current protocol to
enable capturing tunneling in the inverted regime, and addressing the
inaccuracies in rate as we approach classical-limit behavior. While the
current implementation is scalable to atomistic simulations where reactant and
product state are defined in the diabatic framework using methods like the
multi-state empirical valence band model, 84 we are also investigating working
with the adiabatic rather than diabatic states to enable on-the-fly
simulations, an attractive direction given that the MF-RPMD Hamiltonian is
agnostic to the choice of electronic state representation.

##  5 Acknowledgements

N.A., B.A.J., and N. L. acknowledge support from the U.S. Department of
Energy, Office of Basic Energy Sciences, Division of Chemical Sciences,
Geosciences and Biosciences under Award DE-FG02-12ER16362 (Nanoporous
Materials Genome: Methods and Software to Optimize Gas Storage, Separations,
and Catalysis). B.A.J. acknowledges current support by the DOE Office of
Science, Office of Basic Energy Sciences, Division of Chemical Sciences,
Geosciences, and Biosciences, Condensed Phase and Interfacial Molecular
Science program. S.B. acknowledges support from Cornell University, Department
of Chemistry and Chemical Biology.

Additional simulations details for MF-RPMD implementation and details of the
analytic expression for calculating the FGR rate for bilinear coupling model
Hamiltonians are provided in the supporting information.

## References

  * Marcus and Sutin 1985 Marcus, R. A.; Sutin, N. Electron Transfers in Chemistry and Biology. _Biochim. Biophys. Acta_ 1985, _811_ , 265–322
  * Gray and Winkler 1996 Gray, H. B.; Winkler, J. R. ELECTRON TRANSFER IN PROTEINS. _Annu. Rev. Biochem._ 1996, _65_ , 537–561
  * Skourtis et al. 2010 Skourtis, S. S.; Waldeck, D. H.; Beratan, D. N. Fluctuations in Biological and Bioinspired Electron-Transfer Reactions. _Annu. Rev. Phys. Chem._ 2010, _61_ , 461–485
  * Azzouzi et al. 2018 Azzouzi, M.; Yan, J.; Kirchartz, T.; Liu, K.; Wang, J.; Wu, H.; Nelson, J. Nonradiative Energy Losses in Bulk-Heterojunction Organic Photovoltaics. _Phys. Rev. X_ 2018, _8_ , 031055
  * Yarkony 2012 Yarkony, D. R. Nonadiabatic Quantum Chemistry: Past, Present, and Future. _Chem. Rev._ 2012, _112_ , 481–498
  * Hammes-Schiffer and Soudackov 2008 Hammes-Schiffer, S.; Soudackov, A. V. Proton-Coupled Electron Transfer in Solution, Proteins, and Electrochemistry. _J. Phys. Chem. B_ 2008, _112_ , 14108–14123
  * Marcus 1956 Marcus, R. A. On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I. _J. Chem. Phys._ 1956, _24_ , 966–978
  * Wolynes 1987 Wolynes, P. G. Imaginary Time Path Integral Monte Carlo Route to Rate Coefficients for Nonadiabatic Barrier Crossing. _J. Chem. Phys._ 1987, _87_ , 6559–6561
  * Lawrence and Manolopoulos 2018 Lawrence, J. E.; Manolopoulos, D. E. Analytic Continuation of Wolynes Theory into the Marcus Inverted Regime. _J. Chem. Phys._ 2018, _148_ , 102313
  * Fang et al. 2019 Fang, W.; Thapa, M. J.; Richardson, J. O. Nonadiabatic Quantum Transition-State Theory in the Golden-Rule Limit. II. Overcoming the Pitfalls of the Saddle-Point and Semiclassical Approximations. _J. Chem. Phys._ 2019, _151_ , 214101
  * Fang et al. 2020 Fang, W.; Zarotiadis, R. A.; Richardson, J. O. Revisiting Nuclear Tunnelling in the Aqueous Ferrous–Ferric Electron Transfer. _Phys. Chem. Chem. Phys._ 2020, _22_ , 10687–10698
  * Lawrence and Manolopoulos 2020 Lawrence, J. E.; Manolopoulos, D. E. An Improved Path-Integral Method for Golden-Rule Rates. _J. Chem. Phys._ 2020, _153_ , 154113
  * Lawrence and Manolopoulos 2020 Lawrence, J. E.; Manolopoulos, D. E. Confirming the Role of Nuclear Tunneling in Aqueous Ferrous–Ferric Electron Transfer. _J. Chem. Phys._ 2020, _153_ , 154114
  * Lawrence et al. 2019 Lawrence, J. E.; Fletcher, T.; Lindoy, L. P.; Manolopoulos, D. E. On the Calculation of Quantum Mechanical Electron Transfer Rates. _J. Chem. Phys._ 2019, _151_ , 114119
  * Beck et al. 2000 Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H. D. The Multiconfiguration Time-Dependent Hartree (MCTDH) Method: A Highly Efficient Algorithm for Propagating Wavepackets. _Physics Reports_ 2000, _324_ , 1–105
  * Han et al. 2022 Han, S.; Schröder, M.; Gatti, F.; Meyer, H.-D.; Lauvergnat, D.; Yarkony, D. R.; Guo, H. Representation of Diabatic Potential Energy Matrices for Multiconfiguration Time-Dependent Hartree Treatments of High-Dimensional Nonadiabatic Photodissociation Dynamics. _J. Chem. Theory Comput._ 2022, _18_ , 4627–4638
  * Makri 1998 Makri, N. Quantum Dissipative Dynamics: A Numerically Exact Methodology. _J. Phys. Chem. A_ 1998, _102_ , 4414–4427
  * Chatterjee and Makri 2019 Chatterjee, S.; Makri, N. Real-Time Path Integral Methods, Quantum Master Equations, and Classical vs Quantum Memory. _J. Phys. Chem. B_ 2019, _123_ , 10470–10482
  * Makri 2020 Makri, N. Small Matrix Path Integral for System-Bath Dynamics. _J. Chem. Theory Comput._ 2020, _16_ , 4038–4049
  * Tanimura and Kubo 1989 Tanimura, Y.; Kubo, R. Time Evolution of a Quantum System in Contact with a Nearly Gaussian-Markoffian Noise Bath. _J. Phys. Soc. Jpn._ 1989, _58_ , 101–114
  * Tanimura and Wolynes 1991 Tanimura, Y.; Wolynes, P. G. Quantum and Classical Fokker-Planck Equations for a Gaussian-Markovian Noise Bath. _Phys. Rev. A_ 1991, _43_ , 4131–4142
  * Xing et al. 2022 Xing, T.; Li, T.; Liu, Y.; Shi, Q. Mixed Quantum Classical Reaction Rates Based on the Phase Space Formulation of the Hierarchical Equations of Motion. _Chinese Journal of Chemical Physics_ 2022, _35_ , 727–737
  * Redfield 1957 Redfield, A. G. On the Theory of Relaxation Processes. _IBM J. Res. Dev._ 1957, _1_ , 19–31
  * Kimura 2016 Kimura, A. Time-Dependent Renormalized Redfield Theory. _Chemical Physics Letters_ 2016, _645_ , 123–126
  * Tully 1990 Tully, J. C. Molecular Dynamics with Electronic Transitions. _J. Chem. Phys._ 1990, _93_ , 1061–1071
  * Kapral 2006 Kapral, R. Progress in the Theory of Mixed Quantum-Classical Dynamics. _Annu. Rev. Phys. Chem._ 2006, _57_ , 129–157
  * Jain and Subotnik 2015 Jain, A.; Subotnik, J. E. Surface Hopping, Transition State Theory, and Decoherence. II. Thermal Rate Constants and Detailed Balance. _J. Chem. Phys._ 2015, _143_ , 134107
  * Crespo-Otero and Barbatti 2018 Crespo-Otero, R.; Barbatti, M. Recent Advances and Perspectives on Nonadiabatic Mixed Quantum Classical Dynamics. _Chem. Rev._ 2018, _118_ , 7026–7068
  * Miller and Cotton 2016 Miller, W. H.; Cotton, S. J. Classical Molecular Dynamics Simulation of Electronically Non-Adiabatic Processes. _Faraday Discuss._ 2016, _195_ , 9–30
  * Cotton and Miller 2013 Cotton, S. J.; Miller, W. H. Symmetrical Windowing for Quantum States in Quasi-Classical Trajectory Simulations. _J. Phys. Chem. A_ 2013, _117_ , 7190–7194
  * Warburton et al. 2022 Warburton, R. E.; Soudackov, A. V.; Hammes-Schiffer, S. Theoretical Modeling of Electrochemical Proton-Coupled Electron Transfer. _Chem. Rev._ 2022, _122_ , 10599–10650
  * Cao and Voth 1997 Cao, J.; Voth, G. A. A Unified Framework for Quantum Activated Rate Processes. II. The Nonadiabatic Limit. _J. Chem. Phys._ 1997, _106_ , 1769–1779
  * Cao et al. 1995 Cao, J.; Minichino, C.; Voth, G. A. The Computation of Electron Transfer Rates: The Nonadiabatic Instanton Solution. _J. Chem. Phys._ 1995, _103_ , 1391–1399
  * Richardson et al. 2015 Richardson, J. O.; Bauer, R.; Thoss, M. Semiclassical Green’s Functions and an Instanton Formulation of Electron-Transfer Rates in the Nonadiabatic Limit. _J. Chem. Phys._ 2015, _143_ , 134115
  * Richardson 2015 Richardson, J. O. Ring-Polymer Instanton Theory of Electron Transfer in the Nonadiabatic Limit. _J. Chem. Phys._ 2015, _143_ , 134116
  * Lee et al. 2016 Lee, M. K.; Huo, P.; Coker, D. F. Semiclassical Path Integral Dynamics: Photosynthetic Energy Transfer with Realistic Environment Interactions. _Annu. Rev. Phys. Chem._ 2016, _67_ , 639–668
  * Ananth et al. 2007 Ananth, N.; Venkataraman, C.; Miller, W. H. Semiclassical Description of Electronically Nonadiabatic Dynamics via the Initial Value Representation. _J. Chem. Phys._ 2007, _127_ , 084114
  * Church et al. 2018 Church, M. S.; Hele, T. J. H.; Ezra, G. S.; Ananth, N. Nonadiabatic Semiclassical Dynamics in the Mixed Quantum-Classical Initial Value Representation. _J. Chem. Phys._ 2018, _148_ , 102326
  * Liu et al. 2021 Liu, J.; He, X.; Wu, B. Unified Formulation of Phase Space Mapping Approaches for Nonadiabatic Quantum Dynamics. _Acc. Chem. Res._ 2021, _54_ , 4215–4228
  * Cao and Voth 1994 Cao, J.; Voth, G. A. The Formulation of Quantum Statistical Mechanics Based on the Feynman Path Centroid Density. II. Dynamical Properties. _J. Chem. Phys._ 1994, _100_ , 5106–5117
  * Jang and Voth 1999 Jang, S.; Voth, G. A. A Derivation of Centroid Molecular Dynamics and Other Approximate Time Evolution Methods for Path Integral Centroid Variables. _J. Chem. Phys._ 1999, _111_ , 2371–2384
  * Craig and Manolopoulos 2004 Craig, I. R.; Manolopoulos, D. E. Quantum Statistics and Classical Mechanics: Real Time Correlation Functions from Ring Polymer Molecular Dynamics. _J. Chem. Phys._ 2004, _121_ , 3368–3373
  * Craig and Manolopoulos 2005 Craig, I. R.; Manolopoulos, D. E. Chemical Reaction Rates from Ring Polymer Molecular Dynamics. _J. Chem. Phys._ 2005, _122_ , 084106
  * Lawrence and Manolopoulos 2020 Lawrence, J. E.; Manolopoulos, D. E. Path Integral Methods for Reaction Rates in Complex Systems. _Faraday Discuss._ 2020, _221_ , 9–29
  * Novikov et al. 2018 Novikov, I. S.; Suleimanov, Y. V.; Shapeev, A. V. Automated Calculation of Thermal Rate Coefficients Using Ring Polymer Molecular Dynamics and Machine-Learning Interatomic Potentials with Active Learning. _Phys. Chem. Chem. Phys._ 2018, _20_ , 29503–29512
  * Hele and Althorpe 2013 Hele, T. J. H.; Althorpe, S. C. Derivation of a True $(t\rightarrow 0_{+})$ Quantum Transition-State Theory. I. Uniqueness and Equivalence to Ring-Polymer Molecular Dynamics Transition-State-Theory. _J. Chem. Phys._ 2013, _138_ , 084108
  * Althorpe and Hele 2013 Althorpe, S. C.; Hele, T. J. H. Derivation of a True $(t\rightarrow 0_{+})$ Quantum Transition-State Theory. II. Recovery of the Exact Quantum Rate in the Absence of Recrossing. _J. Chem. Phys._ 2013, _139_ , 084115
  * Hele and Althorpe 2016 Hele, T. J. H.; Althorpe, S. C. An Alternative Derivation of Ring-Polymer Molecular Dynamics Transition-State Theory. _J. Chem. Phys._ 2016, _144_ , 174107
  * Habershon et al. 2013 Habershon, S.; Manolopoulos, D. E.; Markland, T. E.; Miller, T. F. Ring-Polymer Molecular Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space. _Annu. Rev. Phys. Chem._ 2013, _64_ , 387–413
  * Menzeleev et al. 2011 Menzeleev, A. R.; Ananth, N.; Miller, T. F. Direct Simulation of Electron Transfer Using Ring Polymer Molecular Dynamics: Comparison with Semiclassical Instanton Theory and Exact Quantum Methods. _J. Chem. Phys._ 2011, _135_ , 074106
  * Kenion and Ananth 2016 Kenion, R. L.; Ananth, N. Direct Simulation of Electron Transfer in the Cobalt Hexammine(II/III) Self-Exchange Reaction. _Phys. Chem. Chem. Phys._ 2016, _18_ , 26117–26124
  * Kretchmer and Miller III 2016 Kretchmer, J. S.; Miller III, T. F. Tipping the Balance between Concerted versus Sequential Proton-Coupled Electron Transfer. _Inorg. Chem._ 2016, _55_ , 1022–1031
  * Liao and Voth 2002 Liao, J.-L.; Voth, G. A. A Centroid Molecular Dynamics Approach for Nonadiabatic Dynamical Processes in Condensed Phases: The Spin-Boson Case. _J. Phys. Chem. B_ 2002, _106_ , 8449–8455
  * Pavese et al. 1999 Pavese, M.; Berard, D. R.; Voth, G. A. Ab Initio Centroid Molecular Dynamics: A Fully Quantum Method for Condensed-Phase Dynamics Simulations. _Chemical Physics Letters_ 1999, _300_ , 93–98
  * Loose et al. 2022 Loose, T. D.; Sahrmann, P. G.; Voth, G. A. Centroid Molecular Dynamics Can Be Greatly Accelerated through Neural Network Learned Centroid Forces Derived from Path Integral Molecular Dynamics. _J. Chem. Theory Comput._ 2022, _18_ , 5856–5863
  * Paesani and Voth 2010 Paesani, F.; Voth, G. A. A Quantitative Assessment of the Accuracy of Centroid Molecular Dynamics for the Calculation of the Infrared Spectrum of Liquid Water. _The Journal of Chemical Physics_ 2010, _132_ , 014105
  * Paesani and Voth 2008 Paesani, F.; Voth, G. A. Nonlinear Quantum Time Correlation Functions from Centroid Molecular Dynamics and the Maximum Entropy Method. _The Journal of Chemical Physics_ 2008, _129_ , 194113
  * Reichman et al. 2000 Reichman, D. R.; Roy, P.-N.; Jang, S.; Voth, G. A. A Feynman Path Centroid Dynamics Approach for the Computation of Time Correlation Functions Involving Nonlinear Operators. _The Journal of Chemical Physics_ 2000, _113_ , 919–929
  * Schwieters and Voth 1998 Schwieters, C. D.; Voth, G. A. The Semiclassical Calculation of Nonadiabatic Tunneling Rates. _J. Chem. Phys._ 1998, _108_ , 1055–1062
  * Hele 2011 Hele, T. J. H. An Electronically Non-Adiabatic Generalization of Ring Polymer Molecular Dynamics. M.Sc. thesis, Oxford University, Oxfordshire, England, 2011
  * Duke and Ananth 2016 Duke, J. R.; Ananth, N. Mean Field Ring Polymer Molecular Dynamics for Electronically Nonadiabatic Reaction Rates. _Faraday Discuss._ 2016, _195_ , 253–268
  * Richardson and Thoss 2013 Richardson, J. O.; Thoss, M. Communication: Nonadiabatic Ring-Polymer Molecular Dynamics. _J. Chem. Phys._ 2013, _139_ , 031102
  * Menzeleev et al. 2014 Menzeleev, A. R.; Bell, F.; Miller, T. F. Kinetically Constrained Ring-Polymer Molecular Dynamics for Non-Adiabatic Chemical Reactions. _J. Chem. Phys._ 2014, _140_ , 064103
  * Kretchmer and Miller III 2016 Kretchmer, J. S.; Miller III, T. F. Kinetically-Constrained Ring-Polymer Molecular Dynamics for Non-Adiabatic Chemistries Involving Solvent and Donor-Acceptor Dynamical Effects. _Faraday Discuss._ 2016, _195_ , 191–214
  * Kretchmer et al. 2018 Kretchmer, J. S.; Boekelheide, N.; Warren, J. J.; Winkler, J. R.; Gray, H. B.; Miller, T. F. Fluctuating Hydrogen-Bond Networks Govern Anomalous Electron Transfer Kinetics in a Blue Copper Protein. _Proc. Natl. Acad. Sci._ 2018, _115_ , 6129–6134
  * Chowdhury and Huo 2017 Chowdhury, S. N.; Huo, P. Coherent State Mapping Ring Polymer Molecular Dynamics for Non-Adiabatic Quantum Propagations. _J. Chem. Phys._ 2017, _147_ , 214109
  * Ananth and Miller 2010 Ananth, N.; Miller, T. F. Exact Quantum Statistics for Electronically Nonadiabatic Systems Using Continuous Path Variables. _J. Chem. Phys._ 2010, _133_ , 234103
  * Ananth 2013 Ananth, N. Mapping Variable Ring Polymer Molecular Dynamics: A Path-Integral Based Method for Nonadiabatic Processes. _J. Chem. Phys._ 2013, _139_ , 124102
  * Duke and Ananth 2015 Duke, J. R.; Ananth, N. Simulating Excited State Dynamics in Systems with Multiple Avoided Crossings Using Mapping Variable Ring Polymer Molecular Dynamics. _J. Phys. Chem. Lett._ 2015, _6_ , 4219–4223
  * Pierre et al. 2017 Pierre, S.; Duke, J. R.; Hele, T. J. H.; Ananth, N. A Mapping Variable Ring Polymer Molecular Dynamics Study of Condensed Phase Proton-Coupled Electron Transfer. _J. Chem. Phys._ 2017, _147_ , 234103
  * Shushkov et al. 2012 Shushkov, P.; Li, R.; Tully, J. C. Ring Polymer Molecular Dynamics with Surface Hopping. _The Journal of Chemical Physics_ 2012, _137_ , 22A549
  * Tao et al. 2018 Tao, X.; Shushkov, P.; Miller, T. F. Path-Integral Isomorphic Hamiltonian for Including Nuclear Quantum Effects in Non-Adiabatic Dynamics. _J. Chem. Phys._ 2018, _148_ , 102327
  * Tao et al. 2019 Tao, X.; Shushkov, P.; Miller, T. F. Simple Flux-Side Formulation of State-Resolved Thermal Reaction Rates for Ring-Polymer Surface Hopping. _J. Phys. Chem. A_ 2019, _123_ , 3013–3020
  * Lawrence and Manolopoulos 2019 Lawrence, J. E.; Manolopoulos, D. E. An Analysis of Isomorphic RPMD in the Golden Rule Limit. _J. Chem. Phys._ 2019, _151_ , 244109
  * Chandler 1987 Chandler, D. _Introduction to Modern Statistical Mechanics_ ; Oxford University Press: New York, 1987
  * Trotter 1959 Trotter, H. F. On the Product of Semi-Groups of Operators. _Proc. Am. Math. Soc._ 1959, _10_ , 545–551
  * Chandler and Wolynes 1981 Chandler, D.; Wolynes, P. G. Exploiting the Isomorphism between Quantum Theory and Classical Statistical Mechanics of Polyatomic Fluids. _The Journal of Chemical Physics_ 1981, _74_ , 4078–4095
  * Craig and Manolopoulos 2005 Craig, I. R.; Manolopoulos, D. E. A Refined Ring Polymer Molecular Dynamics Theory of Chemical Reaction Rates. _J. Chem. Phys._ 2005, _123_
  * Ranya and Ananth 2020 Ranya, S.; Ananth, N. Multistate Ring Polymer Instantons and Nonadiabatic Reaction Rates. _J. Chem. Phys._ 2020, _152_ , 114112
  * Miller et al. 1983 Miller, W. H.; Schwartz, S. D.; Tromp, J. W. Quantum Mechanical Rate Constants for Bimolecular Reactions. _J. Chem. Phys._ 1983, _79_ , 4889–4898
  * Mattiat and Richardson 2018 Mattiat, J.; Richardson, J. O. Effects of Tunnelling and Asymmetry for System-Bath Models of Electron Transfer. _J. Chem. Phys._ 2018, _148_ , 102311
  * Henriksen and Hansen 2008 Henriksen, N. E.; Hansen, F. Y. In _Theories of Molecular Reaction Dynamics: The Microscopic Foundation of Chemical Kinetics_ ; Henriksen, N. E., Hansen, F. Y., Eds.; Oxford University Press, 2008; p 0
  * Lawrence 2020 Lawrence, J. Path-integral methods for non-adiabatic reaction rates. Ph.D. thesis, University of Oxford, Oxfordshire, England, 2020
  * Voth 2006 Voth, G. A. Computer Simulation of Proton Solvation and Transport in Aqueous and Biomolecular Systems. _Acc. Chem. Res._ 2006, _39_ , 143–150