Direct simulation of electron transfer using ring polymer molecular dynamics: Comparison with semiclassical instanton theory and exact quantum methods

The use of ring polymer molecular dynamics (RPMD) for the direct simulation of electron transfer (ET) reaction dynamics is analyzed in the context of Marcus theory, semiclassical instanton theory, and exact quantum dynamics approaches. For both fully atomistic and system-bath representations of condensed-phase ET, we demonstrate that RPMD accurately predicts both ET reaction rates and mechanisms throughout the normal and activationless regimes of the thermodynamic driving force. Analysis of the ensemble of reactive RPMD trajectories reveals the solvent reorganization mechanism for ET that is anticipated in the Marcus rate theory, and the accuracy of the RPMD rate calculation is understood in terms of its exact description of statistical fluctuations and its formal connection to semiclassical instanton theory for deep-tunneling processes. In the inverted regime of the thermodynamic driving force, neither RPMD nor a related formulation of semiclassical instanton theory capture the characteristic turnover in the reaction rate; comparison with exact quantum dynamics simulations reveals that these methods provide inadequate quantization of the real-time electronic-state dynamics in the inverted regime.


I. INTRODUCTION
Condensed-phase electron transfer (ET) reactions are central to many biological and synthetic pathways for energy conversion and catalysis. [1][2][3][4] The development of accurate, robust, and scalable methods for the study of such reactions is thus a key objective in theoretical chemistry. Although transition state theories and rate models for ET have been successfully applied in complex systems, [5][6][7][8] methods for the direct simulation and mechanistic study of ET dynamics in general systems remain less fully developed. To this end, we explore the use of ring polymer molecular dynamics (RPMD) for the description of prototypical ET reactions between mixed-valence transition metal ions in water, and we compare the RPMD approach against benchmark semiclassical and quantum dynamics methods.
Fundamental theoretical challenges in the direct simulation of ET reactions arise due to the coupling of the intrinsically quantum mechanical electronic transitions with slower, classical motions of the surrounding environment. Numerous semiclassical and mixed quantumclassical dynamics methods have been developed for the investigation of electronically non-adiabatic reactions reactions, [9][10][11][12][13][14][15][16][17][18][19] but existing methods do not enable mechanistic studies that are independent of dividing surface assumptions in general systems; nor do they yield dynamical trajectories that preserve the equilibrium Boltzmann distribution 20,21 and allow for the use of rare-event sampling methodologies. 22 New methods are needed to accurately describe coupled electronic and nuclear dynamics and to enable the efficient and robust simulation of long trajectories that bridge the multiple timescales of ET rea) Electronic mail: tfm@caltech.edu actions in complex systems.
RPMD 23 is an approximate quantum dynamical method that is based on the imaginary-time path integral formulation of statistical mechanics. 24,25 It provides an isomorphic classical molecular dynamics model for the real-time evolution of a quantum mechanical system. Previous applications of RPMD include studies of molecular liquids, [26][27][28][29][30][31] hydrogen transfer rates, [32][33][34][35] and tunneling processes in low-dimensional systems. [36][37][38] A key feature of the RPMD method is that it yields realtime molecular dynamics trajectories that preserve the exact quantum Boltzmann distribution and exhibit timereversal symmetry. 23,39 These properties allow RPMD to be used in combination with rare-event sampling methods for the trajectory-based analysis of quantum mechanical tunneling processes in systems involving thousands of atoms. 35,40 We have recently extended RPMD to describe electronic and nuclear dynamics, including solvated electron diffusion 41 and non-adiabatic electron injection into liquid water. 40 In the current paper, RPMD is used to directly simulate ET dynamics in both atomistic and system-bath representations for mixed-valence ET in water. The calculated rates and mechanisms are analyzed in the context of semiclassical and exact quantum methods. A description of the employed methodologies is provided in Section II, and Section III presents the details of the atomistic and system-bath representations. Calculation details are given in Section IV, and a discussion of the results is presented in Section V.

II. METHODS
Several methods are utilized to investigate ET rates and mechanisms, including RPMD, semiclassical instanton theory, exact quantum-mechanical dynamics, and the arXiv:1107.5091v1 [cond-mat.stat-mech] 25 Jul 2011 classical Marcus rate theory for ET. These methods are summarized below.

A. Ring Polymer Molecular Dynamics
The RPMD equations of motion for a quantized electron and N classical particles are 23,41 v (α) = ω 2 n q (α+1) + q (α−1) − 2q (α) where q (α) and v (α) are the position and velocity vectors of α th ring polymer bead, Q j and V j are the position and velocity vectors of the j th classical particle, and n is the number of imaginary-time ring-polymer beads. The intra-bead harmonic frequency is ω n /(β ), where β is the reciprocal temperature. The masses of electron and classical particles are m e and M j , respectively, U ext (q (α) , Q 1 , . . . , Q N ) is the potential energy function of the system, and q (0) = q (n) . Eqs. 1 and 2 generate a classical dynamics that we employ as a model for the real-time dynamics of the system. 41 In the limit of large n, these dynamics preserve the exact Boltzmann distribution. 39 As in classical formulations of the thermal rate constant, [42][43][44] the RPMD rate can be expressed as 32,38 k RPMD = lim t→∞ κ(t)k TST . ( Here, k TST is the transition state theory (TST) approximation for the rate for a dividing surface ξ(r) = ξ ‡ , where ξ(r) is a collective variable, r = q (1) , . . . , q (n) , Q is the full position vector for the system, and Q = {Q 1 , . . . , Q N } denotes the set of classical particle positions. The prefactor, κ(t), is the time-dependent transmission coefficient that accounts for recrossing of trajectories through the dividing surface. An important feature of RPMD is that calculated rates and mechanisms are independent of the choice of TST dividing surface, as in exact quantum and exact classical dynamics. 32,38,45 The TST rate in Eq. 3 is calculated using 33,46,47 Here, F (ξ) is the free energy (FE) along ξ e −β∆F (ξ ) = δ(ξ(r) − ξ ) δ(ξ(r) − ξ r ) , where ξ r is a reference point in the reactant region, and [48][49][50] The scalar r i ∈ {r} in this equation indicates either a ring-polymer or classical particle degree of freedom, m i is the corresponding mass, and d is the total number of degrees of freedom in the system. In Eqs. 4 and 5, . . . denotes the equilibrium ensemble average and . . . c denotes the average in the constrained ensemble Here, where m b is the fictitious Parrinello-Rahman mass, is the full potential energy function for the ring polymer. The transmission coefficient in Eq. 3 is obtained from the flux-side correlation function using where h(ξ) is the Heaviside function,ξ 0 is the initial velocity of the collective variable in an RPMD trajectory released from the dividing surface, and r t is the timeevolved position of the system along that trajectory. 33

B. Semiclassical Instanton Theory
The "Im F" premise in semiclassical rate theory relates the thermal rate constant in the deep-tunneling regime to the analytical continuation of the partition function into the complex plane, [51][52][53][54][55][56] where Q r is the reactant partition function and Im R is the imaginary part of the analytical continuation of the partition function for the full system. In the steepestdescent limit, the "Im F" description is equivalent to the flux-side time correlation formulation 57 of semiclassical instanton (SCI) theory. 58 We adapt this approach to describe the transfer of a single quantized electron in a classical solvent. The partition function for the full system in the ringpolymer representation can be expressed Q n = c dQ I n (Q) (13) where c = N j=1 Mj 2πβ 2 , and is the classical action for a periodic trajectory in imaginary time. The notation presented here assumes that the quantized electron moves in a single dimension. At each solvent configuration, the steepest-descent approximation to I n (Q) is obtained by expanding A({q (α) }; Q) to second order about its global minimum {q (α) }, for which the electron ring-polymer coordinates obey the stationary condition The steepest-descent approximation yields where K is the Hessian matrix given by For a reaction with a barrier, a saddle point satisfies the stationary condition in Eq. 16, and the Hessian matrix exhibits an imaginary normal-mode frequency, η 1 . By analytically continuing η 1 onto the real axis, and by integrating out the zero-frequency normal mode that is associated with the cyclic permutation of the ringpolymer beads, we obtain the steepest-descent SCI rate 58 where (det K) = |η i | 2 is obtained from a product that excludes the zero-frequency mode, and Formal connections between path-integral statistics and reactive tunneling have long been recognized. 25,[58][59][60][61][62] In particular, Althorpe and coworkers 37 have recently emphasized the connection between the TST limit of RPMD and the reversible action work (RAW) formulation of SCI theory. 51,63 To the extent that Eq. 19 is an harmonic approximation to the RAW SCI formulation, 37 where α = 2π(β |η 1 |) −1 , and κ o is the transmission coefficient through a dividing surface that minimizes the recrossing of RPMD trajectories.

C. Exact Quantum Dynamics
We obtain numerically exact quantum dynamics for ET using the Quasi-Adiabatic Path Integral method (QUAPI). [64][65][66][67][68] The method is applied to a redox system composed of two diabatic electronic states and a coordinate representing polarization of the solvent dipole field; the solvent coordinate is in turn linearly coupled to a harmonic oscillator bath.
The Hamiltonian for the redox system 65 where s is the solvent coordinate, p s is the conjugate momentum, and m s is the effective solvent mass. Here, V 11 (s) and V 22 (s) are diabatic states corresponding to reactant and product states for ET, and V 12 (s) is the electronic coupling. The Hamiltonian describing the bath modes and their coupling to the solvent coordinate is where M j , ω j , and Q j are the mass, frequency and the position of the j th bath mode, respectively, and c j is the strength of the coupling between the j th bath mode and the solvent coordinate. The exact quantum mechanical rate constant can be expressed in terms of the symmetrized real-time flux-flux correlation function, 69 where and Q R is the reactant partition function.
Here, H = H S + H B is the full ET Hamiltonian, F = i [H, P 2 ] is the operator for the flux between the reactant and product electronic states, P 2 = |2 2| is the projection operator for the product electronic state, and t c = t − iβ /2 is the complex time. The propagators are discretized into N time slices of length ∆t c , and the trace in Eq. 26 is expanded to yield where Q 0 represents the bath degrees of freedom, s k is the solvent coordinate, and σ k is the electronic state at complex time slice k. The propagators in Eq. 27 are factorized using the quasi-adiabatic short-time approximation 64 Analytical integration over the bath modes then yields where and Here, s 2N +2 = s 1 , s N +2 = s N +1 , and we have introduced the notation s = {s 1 , . . . , s 2N +2 } and σ = {σ 1 , . . . , σ 2N +2 }. In Eq. 31, the path-integral expression for the complextime propagators of the system Hamiltonian is given by The matrix elements in Eq. 32 are obtained using the numerically exact expression where φ m (s, σ) and E m are the eigenstates and eigenenergies of H S , respectively, and M 0 is the number of eigenstates included in the expansion. The discretized form of the non-local influence functional in Eq. 31, which accounts for bath-induced electronic transitions in the system, is where I 0 is the partition function of the uncoupled bath oscillators. 64,70,71 The diagonal elements of {B kk } describe local contributions to the bath response function from a particular complex time slice k along the adiabatic path, and the off-diagonal elements describe nonlocal contributions. For the case of linear system-bath coupling, the diagonal matrix elements are given by and the off-diagonal matrix elements are given by The complex times t k in Eqs. 35 and 36 are provided in Table I.

D. Marcus Theory for ET in a Classical Solvent
In the Marcus theory for ET, 3,72-74 electronic transitions occur at solvent geometries for which the donor and acceptor electronic states are isoenergetic. In the limit of weak electronic coupling and classical solvent motions, the ET rate is thus where V 12 is the electronic coupling matrix element, λ is the solvent reorganization energy, and −∆G 0 is the thermodynamic driving force for the ET reaction. The rate expression in Eq. 37 exhibits three distinct regimes of behavior as the driving force is varied relative to λ. In the normal regime, where −∆G 0 < λ, the rate increases with increasing driving force. A turnover in this trend occurs in the activationless regime, for which −∆G 0 ≈ λ. In the inverted regime, for which −∆G 0 > λ, the rate decreases with increasing driving force.
In the current study, we use implementations for Marcus theory, SCI theory, and RPMD in which the solvent degrees of freedom are treated classically; the role of nuclear quantum effects in diminishing the degree of turnover for the ET rate in the inverted regime 3,75 is not considered here.

III. SYSTEMS
ET dynamics is studied using both all-atom and system-bath representations for mixed-valence transition metal ions in water. These representations are described in this section.

A. Atomistic Representation for ET
The atomistic representation for the ET reaction ( Fig. 1) is described using the potential energy function 76 where q is the electron position and Q is the set of N classical solvent atom positions. Solvent-solvent interactions, U sol (Q), are described using the simple point charge (SPC) model 77 for explicit, rigid water molecules. The remaining interactions are described below, with the values of the parameters provided in Table II.
The electron-water interactions are described using the pairwise pseudopotential 78 where r k = |q − Q k |. For cases in which the atom index k corresponds to a hydrogen atom, and when k corresponds to an oxygen atom, Electron-ion interactions are described using where Q D and Q A denote the respective positions of the donor and acceptor metal ions, which are held fixed at a separation of 6.5Å. These interactions are described using Shaw-type pairwise pseudopotentials. 79 For the acceptor metal ion, where r = |q − Q A |, and for the donor metal ion, with r = |q − Q D |. The asymmetry parameter, , adjusts the thermodynamic driving force for the ET reaction while leaving the solvent reorganization energy unchanged. The values of considered in this study and the corresponding ET regimes are presented in Table III. The ion-water interactions are given by For cases in which atom index k corresponds to a hydrogen atom, and when k corresponds to an oxygen atom, The potential energy functions associated with the acceptor ion, U k A−sol (Q), are obtained by replacing Q D with Q A in Eqs. 46 and 47. These ion-water potential energy functions include electrostatic interactions combined with short-range repulsive terms that reproduce the octahedral coordination structure of the solvated ions. 76

B. System-Bath Representations for ET
The system-bath representation for the ET reaction is described in the position basis using the potential energy function where the first two terms comprise the system potential, and U B is the potential energy contribution due to the bath. The scalar coordinates q and s are the positions of the electron and the solvent mode, respectively.
The first term in the system potential energy function models the ion-electron interaction, (49) This one-dimensional (1D) potential energy function consists of two Coulombic wells capped by parabolic functions to remove the singularity; it is continuous, and its derivative is piecewise continuous over the full range of q. The coefficients in Eq. 49 are provided in Appendix A, and the values of considered for the system-bath representation are presented in Table IV.
The second term in the system potential energy function models the solvent and its interactions with the transferring electron, The first term on the RHS of this equation describes the coupling of the electronic dipole of the redox system to the solvent dipole, and ω s is the effective frequency of the solvent coordinate.
The harmonic oscillator bath potential in Eq. 49 has the same form as in Eq. 24, The bath exhibits Ohmic spectral density with cutoff frequency ω c , where the dimensionless parameter η determines the strength of coupling between the system and the bath modes. 52 The continuous spectral density is discretized into f oscillators with frequencies 32 and coupling constants where j = 1 . . . f . In the current paper, we use two sets of parameters for the system-bath representation. Model SB1 is constructed to reproduce the energy-scales of the atomistic representation, and Model SB2 uses parameters that are numerically less demanding for the QUAPI calculations. The parameters for the models are given in Table V.
As indicated previously, the QUAPI method is implemented using a discrete representation for the diabatic states of the redox system (Eq. 23). The system representation in the position basis described in Eq. 48 is therefore transformed to the electronic diabatic basis for the QUAPI calculations. The resulting diagonal matrix elements for the system potential energy are and and the constant off-diagonal elements V 12 are reported in Table IV. The details of this transformation and the values of the coefficients in Eqs. 55 and 56 are given in Appendix B.

IV. CALCULATION DETAILS
A. Atomistic Representation The atomistic system includes 430 SPC water molecules in a cubic simulation cell with periodic boundary conditions. The side-length of the cell is L = 23.46Å. All calculations are performed at a temperature of T = 300 K, and all pairwise interactions are truncated at a distance of r cut = L/2. Long-range electrostatics are treated by the force-shifting algorithm, 80 where the Coulombic portion of each potential is multiplied by a damping function S(r), such that both the potential and its derivative smoothly vanish at r = r cut . Specifically, Force-shifting reduces the unphysical structuring of water near the the cutoff radius, 80 and it is found to have little effect on the solvent environment of the redox system.

RPMD
The atomistic RPMD simulations are implemented in the DL POLY molecular dynamics package. 81 In all simulations, the RPMD equations of motion are evolved using the velocity Verlet algorithm, 82 and the constraints in the rigid-body water model are implemented using the RATTLE algorithm. 83 The electron is quantized with n = 1024 ring-polymer beads. As in previous RPMD simulations, each timestep for the electron ring polymer involves separate coordinate updates due to forces arising from the physical potential and due to exact evolution of the purely harmonic portion of the ring-polymer potential. The resulting integration algorithm is timereversible and symplectic.
Several collective variables are used to characterize the ET reaction in the atomistic representation. The position of the electron is described by a ring-polymer progress variable, or "bead-count" coordinate, defined as where b = 1.25Å −1 , and the metal ions are symmetrically positioned on the z-axis. We also consider the solvent collective variable where q k ∈ {q H , q O } is the charge on solvent atom k. This solvent collective variable, which is familiar from earlier simulation studies of Marcus theory, 76,84 describes the energy difference between the electonic diabatic states in the tight-binding approximation.
The RPMD rate in Eq. 3 is calculated from the product of the TST rate and the transmission coefficient. The TST rate described in Eq. 4 is obtained from F (f b ), the FE profile in the bead-count coordinate. This FE profile is calculated using umbrella sampling and the weighted histogram analysis method (WHAM), as described below. 85,86 For each value of the asymmetry parameter , the following umbrella sampling protocol is used. Each sampling trajectory is run for at least 50 ps, and thermostatting is performed during the trajectory calculations by resampling the particle velocities from the Maxwell-Boltzmann (MB) distribution every 1.25 ps.
The transmission coefficient in Eq. 11 is calculated using RPMD trajectories that are released from the dividing surface at f ‡ b . For each value of , the dividing surface is chosen to coincide with the maximum along the  Table III. Between 400 and 1200 trajectories are released for each value of . Each RPMD trajectory is evolved for 40 fs with a timestep of 5 × 10 −5 fs and with the initial velocities sampled from the MB distribution. Initial configurations for the released RPMD trajectories are selected every 100 fs from eight long, independent PIMD sampling trajectories that are constrained to the dividing surface ξ ‡ . These sampling trajectories are thermostatted by resampling the velocities every 200 fs, and the constraint to the dividing surface is enforced using the RATTLE algorithm. Two-dimensional (2D) FE surfaces in the ringpolymer centroid coordinate and the solvent coordinate, F (z, ∆U ), are used for the analysis of the ET reaction mechanism. For a given value of , the 2D FE surface is constructed using PIMD sampling trajectories that are harmonically restrained in bothz and ∆U coordinates. Thez coordinate is sampled using 43 uniformly spaced windows in the region of −3.575Å to +3.575Å with a harmonic restraint force constant of 169.7 kcal/molÅ −2 .
To ensure adequate sampling of ring-polymer configurations spanning both metal ions, we use four additional sampling trajectories that are harmonically restrained tō z = ±2.7625Å andz = ±2.925Å with a force constant of 452.5 kcal/molÅ −2 . The solvent coordinate is sampled with 15 uniformly spaced windows in the range −130 to + 150 kcal/mol using a harmonic restraint force constant of 0.023 (kcal/mol) −1 . Each sampling trajectory is run for at least 50 ps, with velocities resampled from the MB distribution every 500 fs. We note that f b is a good progress variable for ET throughout the entire regime of the thermodynamic driving forces, whereas the ring-polymer centroid is not. In the ET inverted regime, the centroid does not fully distinguish between ring-polymer configurations in the reactant and product basins; no such difficulty is experienced in the calculations reported here.

Marcus Theory
Marcus theory rates are calculated using Eqs. 37 and 38. The driving force, −∆G 0 , is obtained from F (∆U ) as the difference between the free energies of the reactant and product minima; these values are reported in Table VI. To the extent that the tight-binding approximation holds, the reorganization energy, λ, is identical for all , and we confirm that this is very nearly the case in our calculations. For the case of symmetric ET ( = 0), the reorganization energy is calculated using λ = 4F (∆U )| ∆U =0 and is found to be 69.7 ± 0.7 kcal/mol.
The coupling matrix element in Eq. 37, |V 12 |, is calculated as 2|V 12 | = E 1 − E 0 , where E 0 and E 1 are the two lowest eigenenergies of the electron in the potential of the isolated metal ions with = 0. These eigenenergies are obtained with an iterative, block Lanczos scheme, 88 performed on a uniform grid of 64 × 64 × 64 points spanning the cubic simulation cell. The iterative Lanczos calculation employs 200 Krylov vectors and an exponential transform parameter of β L = 0.1. The block Lanczos refinement uses ten blocks of five Krylov vectors. This yields a value for the tunnel splitting of |V 12 | = 0.0403 kcal/mol (6.43 × 10 −5 a.u.), which is consistent with previous calculations. 76 This value for the tunnel splitting was assumed to be insensitive to presence of solvent, as has been previously demonstrated, 89 , and independent of the value of the asymmetry parameter . The validity of this latter assumption is confirmed for the system-bath models (see Table IV).

B. System-Bath Representation
As in the atomistic representation, the calculations in the system-bath representation are performed at T = 300 K. The harmonic bath is discretized using f = 12 modes.

RPMD
RPMD rates for the system-bath models are also calculated with the electron quantized using n = 1024 ringpolymer beads. For each value of , the FE profile, F (f b ), is obtained from umbrella sampling along the f b coordinate. For both system-bath models, SB1 and SB2, We note that RPMD results can be affected by coupling of fictitious internal ring-polymer modes to physical frequencies in the system. 90 We thus performed test calculations of the ET rate in these and similar systems using partially adiabatic centroid molecular dynamics (PACMD). 90,91 The PACMD calculations revealed no significant changes from the RPMD results, confirming that this issue does not impact our conclusions.

Marcus Theory
For the calculation of Marcus theory rates, the reorganization energy and the thermodynamic driving force for each value of epsilon are obtained analytically from the diabatic states for the donor and acceptor, V 11 (s) and V 22 (s). For Model SB1, we obtain a solvent reorganization energy of λ = 68.9 kcal/mol, and for Model SB2, we obtain λ = 17.0 kcal/mol. The values of |V 12 | for both system-bath models is given in Table IV.

Semiclassical Instanton Theory
For the SB models, contributions from the linearlycoupled harmonic bath can be factorized and cancelled from the RHS of Eq. 12, yielding expressions that depend only on the electron ring-polymer coordinates and the single classical solvent coordinate, s. Calculation of k SCI then consists of (i) determination of saddle-point configurations for the classical action, A {q (α) }; s , on a numerical grid in the solvent coordinate s, (ii) evaluation of the steepest-descent approximation for I n (s) at each point on the solvent grid, and (iii) integration over the solvent coordinate in Eq. 19 via numerical quadrature. The reactant partition function, Q r , was similarly obtained by evaluating I n (s) via steepest-decent expansion around the minimum-action configuration in the reactant basin. All calculations were performed using n = 2048 beads for the electron ring polymer.
For Model SB1, the grid in the solvent coordinate s consists of 200 uniformly spaced points in the range of −4 to 4 a.u.; for Model SB2, this grid consists of 150 uniformly spaced points in the range −3 to 3 a.u. At each value of s, the saddle-point configuration on the surface A {q (α) }; s corresponds to the maximum along the path of minimum action that connects the reactant and product basins. This path of minimum action is obtained using the string method, 92 with the path discretized into L = 1000 equidistant slices and with minimization performed using Euler integration and a timestep of 2.4 × 10 −3 fs. Initial convergence of the path is achieved when this minimization results in a change of less than 5.3 × 10 −8Å in each degree of freedom. The path is then iteratively refined in the vicinity of the saddle point: a 20-slice sub-section of the path about the saddle point is extracted, the number of slices used to describe the path is doubled, and the sub-section of the path is reminimized with its endpoints fixed. Iterative refinement of the path is complete when the slice of maximum action (i.e., the saddle point configuration) satisfies Eq. 16 to within 10 −5 a.u.

QUAPI
The QUAPI calculation for Model SB2 requires construction of the short-time system propagator followed by two independent Monte Carlo (MC) simulations to evaluate the flux-flux correlation function in Eq. 29.
The complex-time propagator in Eq. 33 is calculated using eigenvalues and eigenfunctions obtained from a 2D discrete variable representation (DVR) grid calculation 93 in the solvent coordinate, s, and the electronic state variable, σ. The DVR Hamiltonian is diagonalized on a grid of 40 uniformly spaced points over a range of −4 to +4 a.u. in s and σ = 1, 2. The number of eigenvalues and eigenvectors used in these calculations (M 0 in Eq. 33) ranges from 30 to 50 for the values of considered in this study.
The flux-flux correlation function in Eq. 29 is obtained from standard path-integral Monte Carlo (PIMC) sampling performed on the 2D DVR grid. In a first PIMC simulation, the correlation function is obtained using where importance sampling is performed using the distribution and the function I i (s, σ; t c ) is defined in Eq. 31. Convergence is achieved with 10 8 MC steps. The normalization constant, D ρ , is obtained from a second, independent PIMC simulation, using Here, importance sampling is performed on the distribution where σ 1 = 2, σ N +1 = 2, σ N +2 = 1, and σ 2N +2 = 1. Convergence is achieved with 10 6 MC steps, and the normalization constant D Λ is obtained by direct matrix multiplication. A maximum of N = 4 path beads are required to converge the flux-flux correlation function over a timescale of 25 fs; no significant changes are observed between calculations performed using N = 4 and N = 8. The reactant partition function is obtained from a single PIMC calculation using the expression where P 1 = |1 1| is the projection operator for the reactant electronic state. The QUAPI calculations for case IV are performed using a larger value for the coupling between the solvent coordinate and the bath modes, η/M ω c = 30. This change  (7) a ET rates are given in s −1 , and the −∆G 0 are given in kcal/mol. The numbers in parentheses denote the statistical uncertainty in the last reported digit.
leads to lower-amplitude oscillations in the flux-flux correlation function and improved numerical convergence of the ET rate calculation. Other features of the fluxflux correlation function, including the timescales for the real-time oscillations and the decorrelation time, are unchanged. The invariance of these features suggests that the parameters used in the current study correspond to the regime in which the ET reaction rate is independent of the solvent-bath coupling. 65,94 RPMD rate calculations performed using different values for η/M ω c also support this conclusion.

A. Atomistic Simulations
The atomistic representation for ET ( Fig. 1) is investigated using direct RPMD simulations and the Marcus rate theory. For each case of the thermodynamic driving force, Fig. 2(a) presents FE profiles for the reactant and product diabatic electronic states as a function of the solvent collective variable, ∆U (Q) (Eq. 59). The FE profiles are obtained by reducing the corresponding 2D surfaces, F (f b , ∆U ), where the reactant and product diabats are associated with ring-polymer configurations for which f b > 0.995 and f b < 0.005, respectively. The results in Fig. 2(a) are graphically identical to those obtained using the tight-binding approximation, and the FE profiles exhibit the anticipated parabolic form, although no assumptions regarding the linear response of the solvent have been made. 76,84 These data, in combination with the calculated tunnel splitting for the transferring electron, are used to calculate the Marcus rates in Table VI. of "kink-pair" configurations, in which the ring polymer spans both redox sites; 25,95,96 a typical kink-pair configuration is illustrated in Fig. 1(a). The dynamical component of the RPMD rate calculation (Eq. 11) is obtained from the long-time plateau 46 of the RPMD transmission coefficient shown in Fig. 2(c). Each transmission coefficient is calculated with respect to a dividing surface at a fixed value of f b , as is described in Sec. IV A. Plateau values in the range of 0.1-0.4 indicate modest recrossing of the RPMD trajectories through these surfaces. For cases in which the thermodynamic driving force corresponds to ET in the normal and acitivationless regimes, Fig. 2(c) illustrates that the RPMD trajectories commit to the reactant or product basins within 10-20 fs, the timescale for local solvent motion between librational rebounds. At thermodynamic driving forces corresponding to the inverted regime, the transmission coefficient plateaus on faster timescales than those involving the rigid solvent molecules. Fig. 3(a) presents a direct comparison of the RPMD and Marcus theory rates throughout the normal and activationless regime for ET in the atomistic representation. The RPMD rates, which are also reported in Table VI, quantitatively agree with the Marcus theory results over 12 orders of magnitude in the ET reaction rate. Unlike the Marcus rates, which are based on a TST description for the reaction, the calculated RPMD rates are independent of any a priori assumptions about the ET reaction mechanism.
Figs. 3(b) and (c) illustrate the ET reaction mechanism that is predicted from the RPMD simulations. Representative RPMD trajectories are projected onto the (z, ∆U ) plane, wherez is the component of the ring-polymer centroid that lies along the axis of the metal ions in the system; also shown are FE profiles for the system in these collective variables. For symmetric ET (Case I), Fig. 3(b) reveals that the RPMD trajectories involve three distinct steps that will be familiar from the Marcus rate theory: (i) solvent fluctuation to a configuration for which the reactant and product diabats are nearly degenerate (indicated by the dashed line), (ii) formation of a kink-pair in the ring-polymer configuration and rapid transfer of the electron from one redox site to the other, and (iii) relaxation of the solvent coordinate in the product basin following the ET event. For ET approaching the activationless regime (Case IV), Fig. 3(c) shows the latter two steps in the mechanism remain, but only a small initial solvent fluctuation is needed to reach solvent configurations for which the electronic diabats are degenerate.
To understand the connection between RPMD and the Marcus theory rate expression, we note that Eq. 37 includes two key terms -an Arrhenius-type contribution that is associated with free energy of solvent reorganization to bring reactant and product diabats into degeneracy and a prefactor that depends on the coupling between the diabatic states. RPMD captures the solvent reorganization energetics because the path-integral-based method preserves exact quantum statistics. 39,87 The RPMD rate and (c) indicate the solvent reorganization mechanism for ET that is anticipated in the Marcus rate theory, and the dashed lines indicate values of ∆U at which the reactant and product diabats cross in Fig. 2(a).
also correctly accounts for the tunneling contribution to the ET reaction rate, which can likewise be attributed to the path-integral basis of the method; the tunnel splitting for the electron between degenerate redox sites is analytically related to the reversible work for forming a kink-pair in the ring-polymer configuration. 25,89,96 Given that the ensemble of reactive RPMD trajectories exhibit the dual rare events of solvent reorganization and kinkpair formation, and given that the FE barriers associated with these two steps are analytically related to the key terms in the Marcus rate expression, it is reasonable that Fig. 3(a) finds good agreement between RMPD and Marcus theory. The RPMD method succeeds in the normal and activationless regimes because it captures the correct physics of the ET reaction. Fig. 4 demonstrates that the success of the RPMD method does not extend into the inverted regime for ET, with both the RPMD rates and reaction mechanism deviating from the predictions of Marcus theory. In Fig. 4(a), the RPMD rates are seen to be only weakly dependent on the increasing driving force, rather than exhibiting the characteristic turnover in this inverted regime. The RPMD trajectories also deviate from the reaction mechanism that is assumed in the Marcus TST, as is seen in Fig. 4(b). The reactive trajectories exhibit kink-pair formation directly from solvent configurations that are characteristic of the reactant basin; the expected solvent reorganization to configurations for which the electronic diabats are degenerate (indicated by the dashed line in the figure) is not observed.
To further explore the successes and failures of RPMD in these various regimes for ET, we compare the method with semiclassical instanton theory and exact quantum dynamics in the following section.

B. System-Bath Simulations
In this section, we employ system-bath representations for ET to allow for the comparison of RPMD with other simulation techniques, including semiclassical instanton and exact quantum dynamics methods. Fig. 5(a) and Table VII present a comparison of the RPMD and Marcus rates for Model SB1, which is parameterized to match the energy-scales for the atomistic representation (Sec. III B). As before, the RPMD method reproduces the Marcus rates throughout the normal and activationless regimes, while failing to predict the turnover of the ET rate in the inverted regime. Analysis of the RPMD reactive trajectories in this system reveals mechanisms that are entirely analogous to those observed in Figs. 3(b), 3(c), and 4(b) for the atomistic system. Specifically, for the normal and activationless regimes, the RPMD trajectories exhibit solvent reorganization to configurations for which the electronic diabats are degenerate, followed by rapid transfer of the electron between redox sites; and for the inverted regime, RPMD predicts ET without prior solvent reorganization. These . The trajectories are plotted as a function of the ring-polymer centroid,z, and the solvent collective variable, ∆U . The FE profile in these collective variables is also presented, with contour lines indicating FE increments of 10 kcal/mol. The white arrows indicate the solvent reorganization mechanism for ET that is anticipated in the Marcus rate theory, and the dashed line indicates the value of ∆U at which the reactant and product diabats cross in Fig. 2(a).
data confirm that Model SB1 exhibits the same essential physics as the atomistic representation. ET rates from the steepest-descent SCI theory (Eq. 19) are also included in Fig. 5(a) and Table VII. Throughout the full range of thermodynamic driving forces, the instanton method tracks the RPMD results, including deviation from the Marcus predictions in the inverted regime. As is shown in Table VII, α-correction of the RMPD rates (Eq. 22, assuming κ o ≈ 1) further improves their agreement with the SCI rates. These results   underscore that the failure of RPMD does not arise from a breakdown in its formal connection with SCI theory; 37 instead, the comparison suggests that both RPMD and the SCI theory share the same underlying flaw in the inverted regime. 97 The mechanistic predictions from SCI theory also show similarities with the RPMD results. Fig. 5(b) presents the Marcus parabolas for the electronic diabats of Model SB1 as a function of the solvent coordinate, s. Also shown are the solvent configurations that correspond to the SCI predictions for the ET transition state. For each value of the thermodynamic driving force, the arrow in the figure indicates the solvent configuration that maximizes I n (s), which corresponds to the largest contribution to the rate in Eq. 19. For the normal and activationless regimes, SCI theory correctly predicts an ET transition state at the crossing of the electronic diabats. However, in the inverted regime, the SCI transition state is instead located at the minimum of the reactant basin. These mechanistic results from SCI theory are consistent with the observed pathways for the RPMD trajectories, which suggests that in the inverted regime, both RPMD and SCI theory overestimate the degree of ET from solvent configurations in the reactant basin.
To further illustrate this issue, we present SCI rate calculations for deep tunneling in a 1D asymmetric double well. Table VIII presents ET reaction rates calculated on the potential energy surface U e-M (q) (Eq. 49), with parameters from Model SB1. Although this is a non-dissipative 1D system, the SCI rate is still welldefined, and it is reported as a function of the potential energy asymmetry. The rates plateau to a finite value with increasing asymmetry, which is consistent with rates for deep tunneling between a bound state and a continuum. [98][99][100] However, this behavior is qualitatively incorrect for tunneling rates between bound states, which should vanish for non-degenerate states in accord with Fermi's Golden Rule. 101 We conclude that SCI theory, as well as the closely related RPMD method, significantly overestimate the tunneling probability between asymmetric bound states, leading to an incorrect ET mechanism and overestimation of the reaction rate in the inverted regime. The ET rates for Model SB1 corresponding to a Marcus-like mechanism (black) and the "direct" mechanism in Eq. 66 (red). SCI rates (blue) correspond to the kinetically favorable mechanism in all regimes. See text for details.
The results for the simple double-well system can be used to deduce a more general argument for the applicability of the RPMD and SCI calculations in ET problems. Table VIII, combined with the condition of detailed balance for the thermal reaction rate, indicates that the SCI rate for transfer in an asymmetric double-well system is approximately For the Marcus-type ET mechanism in which electron tunneling is gated via solvent reorganization that symmetrizes the double-well system, Eq. 65 leads to the TST rate in Eq. 37. However, for an unphysical "direct" ET mechanism in which electron tunneling proceeds from solvent configurations in the reactant basin (i.e., without prior solvent reorganization), Eq. 65 leads to the following TST expression for the ET rate, Fig. 6 presents the ET reaction rates for Model SB1, assuming either the Marcus-type mechanism (Eq. 37, black) or the direct mechanism (Eq. 66, red). Also plotted are the rates calculated using SCI theory (Eq. 19, blue). Throughout the normal and activationless regimes, the rate for the Marcus-type mechanism dominates; in the inverted regime, the rate for the direct mechanism dominates; and the results from SCI theory closely track the larger of these two rates. It is clear that SCI theory (as well as RPMD) features a competition between the correct, Marcus-type mechanism for ET and the unphysical, direct mechanism for ET, and the prevailing mechanism is that which is predicted to be faster. This analysis is fully consistent with the earlier discussions of   (Fig. 5b) for ET in the various regimes. Furthermore, this analysis provides a general basis for expecting the SCI and RPMD methods to accurately describe ET rates in the normal and activationless regimes, and for expecting these methods to significantly overestimate the ET rate in the inverted regime. Table IX presents ET rates for Model SB2, including results obtained using the QUAPI exact quantum dynamics method. Comparison of the RPMD, Marcus theory, and SCI theory rates for ET in Table IX confirms that Model SB2 exhibits all of the previously discussed trends for these approximate methods. Fig. 7 presents the flux-flux correlation functions used to obtain the exact quantum rates for Model SB2.
The results in Fig. 7 emphasize the role of electronic state quantization in the ET reaction dynamics. At larger thermodynamic driving forces, the correlation functions become increasingly oscillatory, with a resonance frequency that matches the electronic state energy gap between the ET reactant and product. 102,103 Integration over this increasingly oscillatory time correlation function (Eq. 25) contributes to the turnover in the ET reaction rate in the inverted regime. The RPMD approximation to the real-time dynamics of the system, which is not expected to capture coherent quantum effects, 23,41 does not fully enforce the quantization of electronic dynamics and leads to the observed inaccuracies in the inverted regime. Approximate quantum dynamical methods that explicitly enforce electronic quantization by using either a discrete electronic state basis or by exactly mapping to a continuous electronic basis are thus expected to provide a better starting point for describing state-to-state electronic dynamics and ET in the inverted regime. 12,14,16,19,104 Further investigation of this point is in progress.

VI. CONCLUSIONS
The current paper demonstrates the applicability of RPMD for the direct simulation of ET reaction dynamics in complex systems. Using both atomistic and systembath representations for ET in a polar solvent, we compare RPMD results with those obtained using Marcus theory, semiclassical instanton theory, and exact quantum dynamics. Throughout the normal and activationless regimes for ET, RPMD correctly predicts the ET reaction mechanism and quantitatively describes the ET reaction rate over 12 orders of magnitude, without invoking any prior mechanistic or transition state assumptions. Analysis of the RPMD trajectories reveals that the accuracy of the method lies in its exact description of statistical fluctuations, with regard to both solvent reorganization and the formation of kink-pair configurations during the electron tunneling event. However, for ET in the inverted regime, both RPMD and SCI theory fail to predict the turnover in the ET reaction rate with increasing thermodynamic driving force. In this regime, both methods overestimate the probability of electronic tunneling from solvent configurations in the reactant basin, leading to an overestimation of the corresponding reaction rates. Exact quantum dynamics calculations illustrate that the limitations of the RPMD method in the inverted regime arise from the inadequate quantization of the real-time electronic-state dynamics; analogous breakdowns of the method have been identified in other applications to strongly coherent quantum systems, including low-dimensional quantum oscillators 23 and electronscattering in dilute fluids. 41 We conclude by emphasizing that the normal and activationless regimes encompass the vast majority ET reactions in biological and synthetic systems. 105 The results presented here thus constitute a significant success for the RPMD method, demonstrating that it allows for the robust, direct simulation of thermally activated ET in systems with over 1000 atoms, leading to the quantitative prediction of ET reaction rates and and the potential discovery and characterization of ET reaction mechanisms in complex systems. A comparable demonstration using other approximate real-time quantum simulation methods has not, to our knowledge, been previously reported. Having established both the applicability and limitations of RPMD for ET reactions dynamics, this work provides the foundation for future studies of ET and proton-coupled ET reactions in enzymes and other condensed-phase systems.

Appendix B: Transformation to Diabatic Basis
The QUAPI method is implemented in an electronic diabatic state representation of the ET reaction. In this appendix, we describe the procedure for transforming the potential energy function for Model SB2 from a position basis for the electron (Eq. 48) to a diabatic basis where the reactant and product electronic states are maximally localized on the donor and acceptor metal atoms.
We begin by calculating the two lowest adiabatic electronic eigenstates (ψ 0 (q; s) and ψ 1 (q; s)) and eigenenergies (E 0 (s) and E 1 (s)) of the system Hamiltonian at fixed values values of the solvent coordinate in the range −8 a 0 ≤ s ≤ 8 a 0 . For each value of s, the system Hamiltonian is diagonalized on a uniform DVR grid of 1024 electron positions in the range −25 a 0 ≤ q ≤ 25 a 0 .
The corresponding potential energy matrix elements in the diabatic basis (Eq. 23) are thus The diagonal elements are found to be parabolic functions of s, and the off-diagonal element are found to be nearly constant with respect to s. We fit V 11 (s) and V 22 (s) to second-order polynomials functions (Eqs. 55 and 56) and employ a constant value for V 12 that corresponds to the s = 0 result. The polynomial expansion coefficients for V 11 (s) and V 22 (s) are provided in Table XIV, and the constant value for V 12 is provided in Table IV.