We investigate the mechanisms of condensed phase proton-coupled electron transfer (PCET) using Mapping-Variable Ring Polymer Molecular Dynamics (MV-RPMD), a recently developed method that employs an ensemble of classical trajectories to simulate nonadiabatic excited state dynamics. Here, we construct a series of system-bath model Hamiltonians for the PCET, where four localized electron-proton states are coupled to a thermal bath via a single solvent mode, and we employ MV-RPMD to simulate state population dynamics. Specifically, for each model, we identify the dominant PCET mechanism, and by comparing against rate theory calculations, we verify that our simulations correctly distinguish between concerted PCET, where the electron and proton transfer together, and sequential PCET, where either the electron or the proton transfers first. This work represents a first application of MV-RPMD to multi-level condensed phase systems; we introduce a modified MV-RPMD expression that is derived using a symmetric rather than asymmetric Trotter discretization scheme and an initialization protocol that uses a recently derived population estimator to constrain trajectories to a dividing surface. We also demonstrate that, as expected, the PCET mechanisms predicted by our simulations are robust to an arbitrary choice of the initial dividing surface.

## I. INTRODUCTION

Understanding the mechanism of condensed phase proton-coupled electron transfer (PCET) in biological processes like water oxidation by photosystem II^{1–4} is an essential step toward the development of biomimetic, photocatalytic materials for water-splitting and the efficient generation of hydrogen fuel.^{5,6} At present, rate theories developed to characterize PCET in specific regimes^{3,7,8} have proven remarkably powerful in explaining experimental observations and more recently in the predicting trends in physical properties and designing catalytic systems.^{8–10} In addition, several direct dynamics simulation methods that can provide mechanistic information have been developed, including on-the-fly coupled electron-nuclear dynamics,^{8,11,12} mixed quantum-classical (MQC) dynamics,^{13–18} and semiclassical simulations.^{19,20} However, these methods employ dynamics that fail to preserve detailed balance and the use of different levels of theory to describe electronic and nuclear motion (particularly by MQC methods) introduce uncontrolled errors in the simulation of nonadiabatic processes.

Imaginary-time path integral^{21} based methods like Ring Polymer Molecular Dynamics (RPMD) overcome this challenge by providing a uniform dynamic framework for electronic and nuclear motion. In addition, RPMD employs an ensemble of classical trajectories that conserve the quantum Boltzmann distribution,^{22} yields reaction rates that are independent of the location of the dividing surface,^{23–26} and can accurately describe the PCET in both adiabatic and nonadiabatic regimes.^{27,28} Unfortunately, RPMD is limited to the simulation of thermal one-electron PCET processes and, in general, cannot be used to characterize nonadiabatic dynamics in multi-level systems.^{29}

Several extensions of RPMD to multi-level systems have been proposed;^{30–34} however, only the Mapping Variable (MV)-RPMD^{31} method incorporates explicit electronic state variables and employs dynamics that conserve the exact quantum Boltzmann distribution. MV-RPMD describes multi-state system dynamics by mapping discrete electronic states to continuous classical analog variables^{35–37} and accurately describes dynamics in both the adiabatic and nonadiabatic regimes.^{31} Recently, two of us demonstrated its short-time accuracy in simulations of photo-induced excited state dynamics in the gas phase.^{34} In this work, we obtain an improved MV-RPMD expression derived from a symmetric rather than asymmetric Trotter discretization scheme,^{38} and we use a recently introduced population estimator^{39} to constrain the ensemble of MV-RPMD trajectories to an arbitrary dividing surface. We then construct a series of system-bath models for the PCET and use MV-RPMD to identify the dominant mechanism in each case. Comparing these mechanistic predictions against rate theory calculations (Fermi’s golden rule for nonadiabatic processes and Kramers rate theory for the adiabatic electron transfer), we show that our simulations correctly distinguish between concerted and sequential PCETs. In addition, we also demonstrate that the mechanistic predictions from MV-RPMD are robust to an arbitrary choice of the dividing surface.

This paper is organized as follows: In Sec. II, we present the modified MV-RPMD expression, and in Sec. III, we describe the quasi-diabatization procedure used to construct system-bath models for the PCET where four localized electron-proton states are coupled to a thermal bath of oscillators via a single solvent coordinate. In Sec. IV, we introduce the MV-RPMD correlation function used to track the electron-proton state population dynamics and the initialization protocol used to constrain MV-RPMD trajectories to a dividing surface. In Sec. V, we provide simulation details, and in Sec. VI, we present the population dynamics obtained from MV-RPMD simulations and we validate the resulting PCET mechanisms against rate theory calculations.

## II. THEORY

### A. MV-RPMD formulation

The Hamiltonian for a general *K*-level system is

where **R** and **P** are nuclear position and momentum operators, respectively, *V*_{0}(**R**) is a state independent nuclear potential, *V*_{nm}(**R**) are elements of the diabatic potential energy matrix, and $\psi n$ represents the *n*th electronic state. Implementing the Meyer-Miller-Stock-Thoss protocol,^{35,36} we map the electronic states to singly excited oscillator (SEO) states,

where $an\u2020$ and *a*_{m} are boson creation and annihilation operators, respectively, that obey the commutation rules $[an\u2020,am]=\delta nm$. In Eq. (2), we use the notation $n=0102\u20261n\u20260K$, to represent SEO states that correspond to a product of *K* − 1 uncoupled oscillators in the ground state and one oscillator in the first excited state.

Following the original MV-RPMD derivation,^{31} path integral discretization of the canonical partition function, $Z=Tre\u2212\beta \u0124$, where *β* = 1/*kT*, is performed using continuous Cartesian variables for the electronic and nuclear degrees of freedom by inserting *N* − 1 copies of the identity,^{37}

where $P\u2261\u2211nnn$ is the projection operator in the SEO basis. Evaluating the matrix elements of the Boltzmann operator using the symmetric Trotter approximation (detailed derivation provided in Appendix A) and employing a Wigner transform in the electronic variables,^{31} we obtain an exact path integral expression for the quantum Boltzmann distribution in electronic and nuclear phase space variables,

where we set $\u210f=1$ throughout the manuscript, *β*_{N} = *β*/*N*, $\u222bd{R\alpha}\u2009\u2261\u2009\u222bdR1\u222bdR2\u2026\u222bdRN$, and similarly for the other variables of integration. In Eq. (4), the MV-RPMD Hamiltonian is

where *N* is the number of ring polymer beads, and the nuclear ring polymer Hamiltonian is

where *M* is the physical mass of the nuclei and *ω*_{N} = *N*/*β*. The electron-nuclear interaction term in Eq. (5) is

where

and **x**_{α} and **p**_{α} are continuous position and momentum vectors of length *K* representing the *K* electronic states of the *α*th ring polymer bead. Finally, the interaction matrix in Eq. (8) is given by

a result that is well known in the context of state space path integrals.^{40} The interaction matrix in Eq. (10) is symmetric (in keeping with the original quantum Hamiltonian) making the MV-RPMD Hamiltonian symmetric and improving the numerical stability of the approximate dynamics. We also emphasize that the symmetric and asymmetric formulations are equivalent for equilibrium simulations and exhibit similar bead-convergence properties.

## III. PCET MODEL SYSTEMS

Previous work using RPMD for the simulation of the PCET in condensed phase model systems used a position-space representation to describe a single distinguishable electron and proton coupled to a thermal bath.^{27} Exact quantum dynamics studies^{41} and surface hopping based simulations^{13} for similar model systems choose to employ a two-state representation of the electron donor and acceptor states coupled to a position space proton. Here, we transform these model Hamiltonians to a representation where four localized, quasi-diabatic electron-proton states are coupled to a thermal bath via a solvent polarization coordinate. The quasi-diabatic states are labeled, DD, DA, AD, and AA following previous literature,^{13} where the letters D/A indicate the donor/acceptor state of the particle and the first letter describes the state of the electron, while the second letter describes the state of the proton.

Following the quasi-diabatization procedure (described in detail in Appendix B), we obtain a four-state system-bath PCET Hamiltonian,

where *s*, *P*_{s}, and *m*_{s} are the position, momentum, and mass of the solvent polarization coordinate and $VXY,X\u2032Y\u2032(s)$ are the elements of the diabatic potential energy matrix where the subscripts $X/Y/X\u2032/Y\u2032={D,A}$ label the donor and acceptor states of the particles. In Eq. (11), *P*_{j}, *Q*_{j}, and *M* are the momentum, position, and mass of the *j*th bath mode, and *c*_{j} is the coupling between the solvent and the *j*th bath mode of frequency *ω*_{j}. The bath spectral density is Ohmic,

with cut-off frequency *ω*_{c} = *ω*_{s}, and the dimensionless parameter *η*/*m*_{s}*ω*_{s} determines the coupling strength between the solvent and the bath modes.^{42} The continuous spectral density is discretized into *f* oscillators with frequencies^{23}

and the coupling constants *c*_{j} are defined as

where *j* = 1, *…*, *f*.

The diagonal elements of the potential energy matrix in Eq. (11) obtained through our quasi-diabatization protocol are fitted to quadratic polynomials of the form

and the off-diagonal couplings are taken to be constants that are independent of the solvent coordinate.

## IV. STATE POPULATION DYNAMICS

In general, thermal real-time correlation functions in the MV-RPMD framework are written as

where ${\xi \alpha}t$ represents the set of bead positions and momenta {**R**_{α}, **P**_{α}, **x**_{α}, **p**_{α}} at time *t*, and the bead-averaged functions $A({\xi \alpha}0)=1/N\u2211\alpha A(\xi \alpha (0))$ and $B({\xi \alpha}t)$ is similarly defined. The initial positions and momenta are generated from a standard Path Integral Monte Carlo (PIMC) simulation that employs the sampling function *W*. For a system initially at equilibrium, $W=e\u2212\beta NHN({\xi \alpha}0)$ with the MV-RPMD Hamiltonian, *H*_{N}, defined in Eq. (5); however, this function can also be defined to describe an initial non-equilibrium distribution as discussed below. Real-time trajectories are generated by integrating equations of motion corresponding to the MV-RPMD Hamiltonian,

For the PCET model systems considered here, the nuclear position vector, **R**_{α} = (*s*_{α}, **Q**_{α}), includes both the 1D solvent coordinate coupled to the local electron-proton states and the positions of all the bath modes.

Here, we investigate the mechanism of thermal PCET by initializing trajectories to a non-equilibrium distribution, *ρ*_{neq}(0), corresponding to a particular choice of the dividing surface. We then track the electron-proton state population dynamics by evaluating the real-time quantum correlation function,

where the Heaviside function, *h*, is defined in terms of the solvent coordinate and allows us to separately ensemble average over trajectories moving forward (from the dividing surface toward reactants) and backwards (toward products),

In the MV-RPMD framework, the Heaviside function in Eq. (18) is written in terms of the solvent ring polymer centroid, $h\u2261h(\xb1(s\xaft\u2212s\u2021))$, where $s\xaf=1/N\u2211\alpha =1Ns\alpha $. The *n*th state populations at time *t* are evaluated using the “Boltzmann” estimator,^{31,34}

where **Γ**_{nn} is a diagonal element of the matrix previously defined in Eq. (8) and the time-evolved positions and momenta are obtained by integrating the MV-RPMD equations of motion in Eq. (17).

To initialize trajectories to the dividing surface, we define an initial non-equilibrium density operator, $\rho neq\u2009=\u2009\rho neqsys\u2009\u2297\u2009\rho eqbath$, where the full system is divided into a relevant subsystem described with non-equilibrium initial conditions and the bath that is initially at equilibrium. The subsystem density matrix is defined by

where *H*_{s} is the subsystem Hamiltonian given by the first line of Eq. (11), *P*_{n} is the population of the *n*th state, and the solvent position, $s\u2021$, and electron-proton state populations, $Pn\u2021$, together define the dividing surface. Ignoring the Boltzmann weights associated with each electronic state, we can write the corresponding constraints in the MV-RPMD framework as

where the nuclear ring polymer Hamiltonian is defined in Eq. (6) and $s\xaf0$ is the nuclear RP centroid constrained to its dividing surface value, $s\u2021$. Further, in Eq. (22), we use the recently derived “semiclassical” estimator,^{39}

where $PnSC\alpha $ is the state population associated with the *α*th bead. We note that this population estimator was rigorously derived to yield the exact equilibrium populations at time *t* = 0^{39} and is of similar form to the original semiclassical population function.^{35,36} The present bead-averaged form in Eq. (23) has also been used as an estimator in the nonadiabatic-RPMD method where trajectories are initialized to an exact equilibrium path-integral distribution and time-evolved under the semiclassical mapping Hamiltonian.^{43} Finally, it is important to recognize that constraining electronic state populations via $PnSC$ in the correlation function in Eq. (22) does not constrain $Pn\beta $ to the same values at *t* = 0 since the latter includes the correct Boltzmann weights for each electronic state at a given nuclear configuration.

## V. SIMULATION DETAILS

We construct three model systems that correspond to different PCET regimes and report values of shared parameters for each case in Table I. Parameters for the quasi-diabatic potential energy matrix elements are tabulated in Appendix C.

Parameter^{a}
. | Model I . | Model II . | Model III . |
---|---|---|---|

m_{s} | 22 000 | 22 000 | 22 000 |

ω_{s} × 10^{4} | 3.72 | 4.00 | 3.72 |

f | 12 | 12 | 12 |

M | m_{s} | m_{s} | m_{s} |

η/m_{s}ω_{s} | 1 | 1 | 1 |

T/K | 300 | 300 | 300 |

Parameter^{a}
. | Model I . | Model II . | Model III . |
---|---|---|---|

m_{s} | 22 000 | 22 000 | 22 000 |

ω_{s} × 10^{4} | 3.72 | 4.00 | 3.72 |

f | 12 | 12 | 12 |

M | m_{s} | m_{s} | m_{s} |

η/m_{s}ω_{s} | 1 | 1 | 1 |

T/K | 300 | 300 | 300 |

^{a}

All parameters specified in atomic units.

For each model, we calculate the real-time correlation function in Eq. (18) by sampling the initial nuclear and electronic non-equilibrium distribution using Path Integral Monte Carlo (PIMC). The initial electronic state population variables should be sampled subject to the bead-average constraint described in Eq. (22). However, following previous work,^{34} we implement this constraint by setting *individual* bead state populations to the desired values at the dividing surface rather than constraining the average,

The dividing surface for all three models is chosen to be the intersection of the reactant (DD) and product (AA) quasi-diabatic state potentials such that $s\u2021$ = 0 a.u., and only the DD and AA states are populated with $PDD\u2021\u2009=\u2009PAA\u2021\u2009=\u20090.5$ and $PDA\u2021=\u2009PAD\u2021\u2009=\u20090$. For each model, we sample the distribution with a total of 5 × 10^{8} MC points and bead convergence is achieved with *N* = 10 beads.

For all three models, MV-RPMD trajectories initialized to the dividing surface are propagated using a 4th order Adams-Bashforth-Moulton predictor corrector integrator with a time step of size 2.4 × 10^{−2} fs. Trajectories were integrated for a total simulation time of 1200 fs for models I and III and 9600 fs for model II. The number of trajectories used to obtain the converged results shown below was 2.5 × 10^{4}, 8 × 10^{4}, and 1.5 × 10^{5} for models I, II, and III, respectively.

We separate the ensemble of trajectories into a group that moves “forward” toward product formation (increasing values of the solvent coordinate) and a group that moves “backward” toward the reactant state (decreasing values of the solvent coordinate) to obtain the correlation function $CPn,h(t)$ defined in Eq. (18). Splicing the forward and backward averages together at time zero, we obtain the population plots shown here.

Finally, we use model III to demonstrate that the mechanism predicted by MV-RPMD is independent of the choice of the initial dividing surface. We choose a different dividing surface with $s\u2021$ = −0.8 a.u. (at the intersection of the DD and AD states), and the initial electronic state populations are taken to be $PDD\u2021\u2009=\u2009PAD\u2021\u2009=\u20090.5$ and $PDA\u2021\u2009=\u2009PAA\u2021\u2009=\u20090$. For this simulation, trajectories were integrated for a total time of 1200 fs and 2.5 × 10^{4} trajectories were employed to obtain the converged results shown here.

## VI. RESULTS AND DISCUSSION

The diabatic potential energy surfaces as a function of the solvent coordinate for model I are shown in Fig. 1, and the corresponding population dynamics are shown in Fig. 2. Reading the plot chronologically from left to right, we find the initially populated reactant state (DD) where both electron and proton are in the donor state transfers population to the product state (AA) where both the electron and proton are in the acceptor state. This indicates a concerted PCET mechanism where the proton and electron transfer simultaneously on a sub-picosecond time scale. The energetically unfavorable AD and DA states are not involved in the PCET process, but we find a small population in both states that decays to zero at long times.

We plot the diabatic potential energy surfaces for model II in Fig. 3, and the corresponding MV-RPMD population dynamics are plotted in Fig. 4. Again, reading the plot chronologically, we find that both the reactant (DD) state and the DA (proton transfer only) state are populated at long times, *t* = −3000 fs.

In Fig. 4, we see additional population transfer from the DD to DA state on a timescale of ≈4800 fs preceding the rise in the product (AA) state population. We also note a negligible population transfer from the DA to AD state at short times that decays into thermal population in the AD state at longer times. These results thus suggest a sequential mechanism for the PCET where the proton transfers first, facilitating the electron transfer.

The diabatic potential energy surfaces for model III is shown in Fig. 5, and the corresponding population dynamics are shown in Fig. 6. We find that the system is initially in the reactant DD state with significant thermal population in the DA state. Following the dynamics, we find, however, that unlike model II population transfers from the reactant state to the AD state corresponding to the electron transfer preceding the rise in population of the product AA state. This indicates a sequential PCET mechanism where the electron transfers first facilitating the proton transfer.

Despite initializing MV-RPMD trajectories to the same initial dividing surface for all three models, we find population dynamics point to three different PCET mechanisms. We now show that MV-RPMD simulations can yield mechanistic insights independent of the initial choice of the dividing surface for the reactive trajectories by using a different dividing surface in model III. In Fig. 7, we plot the results of this simulation where the initial dividing surface is chosen to be at the intersection of the reactant DD state and the electron-transfer only at the AD state. We find that the predicted mechanism is unchanged—population transfer from the reactant state to the AD state first, before PCET product formation.

### A. Verification with rate theories

We verify the accuracy of the PCET mechanism predicted by the MV-RPMD simulation by calculating Fermi’s Golden Rule (FGR) rates for the concerted PCET, electron-transfer, and proton-transfer for each model.^{44} For Models I and III, the electron transfer is near-adiabatic and we use Kramers rate theory^{45} to calculate rates for these processes.

We estimate the FGR rate using a simple analytical form derived for systems in which the reactant and product diabatic potential energy surfaces are displaced harmonic oscillators with frequency *ω* and coupling Δ,^{46,47}

where $\omega =2a/ms$ is the frequency of the product diabatic state, $v$ = (*V*_{R} − *V*_{P})/*ω*, *z* = *βω*/2, $S=ms\omega Vd2/2\u210f$, $Iv$ is a modified Bessel function of the first kind, *V*_{d} is the horizontal displacement of the diabatic potential energy functions, and *V*_{R/P} are the values of the potential energy at the reactant/product minimum such that *V*_{R} − *V*_{P} measures the driving force. For adiabatic ET, we use Kramers theory,^{45}

where *ω*_{b} is the frequency at the top of the barrier, $Gcl\u2021$ is the solvent FE barrier when the solvent is treated classically, and *γ* = *η*/*M*_{S}.^{48} The resulting rates are reported in Table II, and as expected, we find that the fastest rate for model I corresponds to a concerted PCET reaction, for model II, the proton transfer reaction is the most rapid, and for model III, the electron transfer reaction rate is the fastest.

Reaction path . | Model I . | Model II . | Model III . |
---|---|---|---|

k_{DD→AA} | 1.85 × 10 ^{7} | 1.61 × 10^{6} | 4.70 × 10^{6} |

k_{DD→DA} | 9.81 × 10^{−17} | 2.53 × 10 ^{9} | 5.97 × 10^{4} |

k_{DD→AD} | 2.69 × 10^{5*} | 1.01 × 10^{6} | 1.03 × 10 ^{11*} |

Reaction path . | Model I . | Model II . | Model III . |
---|---|---|---|

k_{DD→AA} | 1.85 × 10 ^{7} | 1.61 × 10^{6} | 4.70 × 10^{6} |

k_{DD→DA} | 9.81 × 10^{−17} | 2.53 × 10 ^{9} | 5.97 × 10^{4} |

k_{DD→AD} | 2.69 × 10^{5*} | 1.01 × 10^{6} | 1.03 × 10 ^{11*} |

## VII. CONCLUSIONS

We have extended the applicability of MV-RPMD to the simulation of the condensed phase PCET using an improved formalism and a new population estimator to follow state-to-state population transfer dynamics. We employed a simple quasi-diabatization procedure to build three model PCET systems where four local electron-proton states are coupled to a thermal bath via a single solvent polarization coordinate. Following the population dynamics by initializing MV-RPMD trajectories to an arbitrary dividing surface, we identify the mechanism of the PCET for each of the three models and verify the accuracy of the predicted mechanism against FGR and Kramers rate theory predictions. By performing a simulation with a different dividing surface, we were also able to clearly establish that our MV-RPMD simulations yield mechanisms that are independent of the initial choice of the dividing surface to which trajectories are constrained.

The direct dynamic simulation techniques presented here can be readily extended to future studies of complex photochemical reactions and particularly photo-initiated PCET processes in the condensed phase. Future work in this direction will include deriving a systematic correction to the approximate MV-RPMD dynamics. In addition, we recognize that accurately parameterizing a system-bath Hamiltonian of the form described in Appendix B from an atomistic simulation remains a significant challenge.

## ACKNOWLEDGMENTS

The authors acknowledge funding support for this work from the National Science Foundation CAREER Award No. CHE-1555205. N.A. acknowledges partial funding support from an Alfred P. Sloan Foundation Fellowship and a Cottrell Scholar Award. S.P. acknowledges support from an NDSEG Fellowship No. FA9550-11-C-0028 awarded by the Department of Defense, Air Force Office of Scientific Research. J.R.D. acknowledges funding support from a National Science Foundation Graduate Research Fellowship under Grant No. DGE-1144153.

### APPENDIX A: SYMMETRIC TROTTER DERIVATION

In the limit that *N* → *∞*, the high-temperature symmetric Trotter approximation is used to separate the state independent nuclear potential operator, *V*_{0}, and the diabatic potential energy matrix, *V*, from the nuclear kinetic energy operator *T*,

The nuclear kinetic matrix element can be evaluated exactly to obtain

Substituting Eq. (A2) back in the Boltzmann matrix element, we have

In order to evaluate the electronic matrix element, we begin by defining a diagonal matrix with elements, $VDR\alpha ,R\alpha +1\u2009=\u200912(VD(R\alpha )\u2009+\u2009VD(R\alpha +1))$, and off-diagonal matrix elements $VODR\alpha ,R\alpha +1\u2009=\u200912(VOD(R\alpha )+VOD(R\alpha +1))$. Employing a high-temperature Trotter approximation, we further split the off-diagonal terms symmetrically around the diagonal terms to obtain

The off-diagonal matrix elements are easily evaluated,

where *V*_{nj} is used to indicate off-diagonal elements of the diabatic potential energy matrix. Substituting Eq. (A5) into Eq. (A4), we obtain an expression for the electronic matrix elements by considering two cases:

(

*n*=*m*)**:**If

*n*=*j*,

If

*n*$\u2260$*j*,

(

*n*$\u2260$*m*)**:**If

*n*=*j*and if*m*$\u2260$*j*,

If

*n*$\u2260$*j*and if*m*=*j*,

If

*n*$\u2260$*j*and if*m*$\u2260$*j*,

### APPENDIX B: QUASI-DIABATIZATION PROTOCOL

To construct a Hamiltonian in the basis of local electron-proton states, we start with a previously-used system-bath model Hamiltonian for the PCET where the proton is represented in position space and a two-state system describes the electron transfer.^{13,41,49–51} The system Hamiltonian is

In Eq. (B1), *R* is the proton coordinate with conjugate momentum *P*_{R}, and *V*_{p}(*R*) is a double well potential in the proton coordinate,

where *m*_{R} is the mass of the proton, *ω*_{R} is the frequency, *λ* is a measure of anharmonicity, and *V*_{0} determines the height of the barrier for the proton transfer. Further, the proton-solvent coupling is

where *μ*_{1} and *ϕ* are constants that can be chosen to favor either concerted or sequential mechanism. The two-state diabatic potential for the electron transfer is

where *μ*_{2}, *a*_{i}, and *ϕ* are constants that can be tuned to construct models that favor either concerted or sequential mechanisms. Parameters for the three models considered here are provided in Table III.

Parameter^{a}
. | Model I . | Model II . | Model III . |
---|---|---|---|

m_{R} | 1836.1 | 1836.1 | 1836.1 |

ω_{R} | 0.0104 | 0.0104 | 0.0104 |

V_{0} | 0.012 | 0.014 | 0.012 |

s_{1} | −2.13 | −2.16 | −2.13 |

s_{2} | 2.13 | 2.16 | 2.13 |

V_{12} | 0.002 45 | 0.0124 | 0.002 45 |

μ_{1} | 0.001 1 | 0.017 | −0.001 1 |

μ_{2} × 10^{3} | 5.84 | 0.71 | 5.84 |

λ | 0.0 | 0.012 | 0.0 |

Parameter^{a}
. | Model I . | Model II . | Model III . |
---|---|---|---|

m_{R} | 1836.1 | 1836.1 | 1836.1 |

ω_{R} | 0.0104 | 0.0104 | 0.0104 |

V_{0} | 0.012 | 0.014 | 0.012 |

s_{1} | −2.13 | −2.16 | −2.13 |

s_{2} | 2.13 | 2.16 | 2.13 |

V_{12} | 0.002 45 | 0.0124 | 0.002 45 |

μ_{1} | 0.001 1 | 0.017 | −0.001 1 |

μ_{2} × 10^{3} | 5.84 | 0.71 | 5.84 |

λ | 0.0 | 0.012 | 0.0 |

^{a}

All parameters specified in atomic units.

For each value of the solvent configuration in the range −6*a*_{0} ≤ *s* ≤ 6*a*_{0}, we diagonalize the system Hamiltonian on a uniform Discrete Variable Representation (DVR) grid in the proton coordinate with a grid range of −2*a*_{0} ≤ *R* ≤ 2*a*_{0} and 100 grid points. The adiabatic eigenstates obtained upon diagonalizing the system Hamiltonian are written as $R;s|\u03f5i$, where *ϵ*_{i} is the *i*th adiabatic state with eigenenergy *E*_{i}.

Further, by diagonalizing the system Hamiltonian for a single electronic state (donor or acceptor) at each value of *s*, we construct localized proton wavefunctions, $\u27e8R;s|lj\u27e9$, where *l*_{j} is the *j*th quasi-diabatic local electron-proton states that can be expressed in terms of the adiabatic eigenstates as

For a given *s* value, matrix elements of the Hamiltonian in the quasi-diabatic basis can then be constructed using

where *E*_{i} is the energy of the *i*th eigenstate of the Hamiltonian in Eq. (B1).

The overlap between the reference quasi-diabatic wavefunction and the adiabatic state for a given value of the solvent coordinate, *s*, is then obtained by evaluating

### APPENDIX C: PARAMETERS FOR QUASI-DIABATIC POTENTIAL SURFACES

We provide the diabatic potential energy matrix parameters for all three models below (Tables IV–IX).

Diabatic . | a
. | b
. | c
. |
---|---|---|---|

V_{DD} | 0.0015 | 0.0075 | −0.0041 |

V_{DA} | 0.0015 | 0.0055 | 0.0072 |

V_{AD} | 0.0015 | −0.0055 | 0.0072 |

V_{AA} | 0.0015 | −0.0075 | −0.0041 |

Diabatic . | a
. | b
. | c
. |
---|---|---|---|

V_{DD} | 0.0015 | 0.0075 | −0.0041 |

V_{DA} | 0.0015 | 0.0055 | 0.0072 |

V_{AD} | 0.0015 | −0.0055 | 0.0072 |

V_{AA} | 0.0015 | −0.0075 | −0.0041 |

Coupling . | Δ . |
---|---|

V_{DD,DA} | 9.7 × 10^{−5} |

V_{DD,AD} | 2.5 × 10^{−3} |

V_{DD,AA} | 1.8 × 10^{−4} |

V_{DA,AD} | 1.8 × 10^{−4} |

V_{DA,AA} | 2.5 × 10^{−3} |

V_{AD,AA} | 9.7 × 10^{−5} |

Coupling . | Δ . |
---|---|

V_{DD,DA} | 9.7 × 10^{−5} |

V_{DD,AD} | 2.5 × 10^{−3} |

V_{DD,AA} | 1.8 × 10^{−4} |

V_{DA,AD} | 1.8 × 10^{−4} |

V_{DA,AA} | 2.5 × 10^{−3} |

V_{AD,AA} | 9.7 × 10^{−5} |

Diabatic . | a
. | b
. | c
. |
---|---|---|---|

V_{DD} | 0.0015 | 0.0072 | −0.0018 |

V_{DA} | 0.0018 | 0.0058 | −0.0013 |

V_{AD} | 0.0018 | −0.0061 | 0.0034 |

V_{AA} | 0.0016 | −0.0083 | −0.0018 |

Diabatic . | a
. | b
. | c
. |
---|---|---|---|

V_{DD} | 0.0015 | 0.0072 | −0.0018 |

V_{DA} | 0.0018 | 0.0058 | −0.0013 |

V_{AD} | 0.0018 | −0.0061 | 0.0034 |

V_{AA} | 0.0016 | −0.0083 | −0.0018 |

Coupling . | Δ . |
---|---|

V_{DD,DA} | 1.1 × 10^{−3} |

V_{DD,AD} | 1.2 × 10^{−4} |

V_{DD,AA} | 1.2 × 10^{−4} |

V_{DA,AD} | 1.2 × 10^{−4} |

V_{DA,AA} | 1.2 × 10^{−4} |

V_{AD,AA} | 1.4 × 10^{−3} |

Coupling . | Δ . |
---|---|

V_{DD,DA} | 1.1 × 10^{−3} |

V_{DD,AD} | 1.2 × 10^{−4} |

V_{DD,AA} | 1.2 × 10^{−4} |

V_{DA,AD} | 1.2 × 10^{−4} |

V_{DA,AA} | 1.2 × 10^{−4} |

V_{AD,AA} | 1.4 × 10^{−3} |

Diabatic . | a
. | b
. | c
. |
---|---|---|---|

V_{DD} | 0.0015 | 0.008 | 0.0009 |

V_{DA} | 0.0015 | 0.0098 | 0.013 |

V_{AD} | 0.0015 | −0.0056 | −0.0095 |

V_{AA} | 0.0015 | −0.013 | 0.0009 |

Diabatic . | a
. | b
. | c
. |
---|---|---|---|

V_{DD} | 0.0015 | 0.008 | 0.0009 |

V_{DA} | 0.0015 | 0.0098 | 0.013 |

V_{AD} | 0.0015 | −0.0056 | −0.0095 |

V_{AA} | 0.0015 | −0.013 | 0.0009 |

Coupling . | Δ . |
---|---|

V_{DD,DA} | 6.9 × 10^{−4} |

V_{DD,AD} | 2.5 × 10^{−3} |

V_{DD,AA} | 1.8 × 10^{−4} |

V_{DA,AD} | 1.8 × 10^{−4} |

V_{DA,AA} | 2.5 × 10^{−3} |

V_{AD,AA} | 6.9 × 10^{−4} |

Coupling . | Δ . |
---|---|

V_{DD,DA} | 6.9 × 10^{−4} |

V_{DD,AD} | 2.5 × 10^{−3} |

V_{DD,AA} | 1.8 × 10^{−4} |

V_{DA,AD} | 1.8 × 10^{−4} |

V_{DA,AA} | 2.5 × 10^{−3} |

V_{AD,AA} | 6.9 × 10^{−4} |