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.

Understanding the mechanism of condensed phase proton-coupled electron transfer (PCET) in biological processes like water oxidation by photosystem II1–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 regimes3,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 integral21 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)-RPMD31 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 variables35–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 estimator39 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.

The Hamiltonian for a general K-level system is

Ĥ=PTP2M+V0(R)+n,m=1KψnVnm(R)ψm,
(1)

where R and P are nuclear position and momentum operators, respectively, V0(R) is a state independent nuclear potential, Vnm(R) are elements of the diabatic potential energy matrix, and ψn represents the nth electronic state. Implementing the Meyer-Miller-Stock-Thoss protocol,35,36 we map the electronic states to singly excited oscillator (SEO) states,

|ψnψm|anam|nm|,
(2)

where an and am are boson creation and annihilation operators, respectively, that obey the commutation rules [an,am]=δnm. In Eq. (2), we use the notation n=01021n0K, 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βĤ, 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 

I=dxdR|x,Rx,R|P,
(3)

where Pnnn 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,

ZlimNd{Rα}d{Pα}d{xα}d{pα}×eβNHN({Rα},{Pα},{xα},{pα})sgn(Θ),
(4)

where we set =1 throughout the manuscript, βN = β/N, d{Rα}dR1dR2dRN, and similarly for the other variables of integration. In Eq. (4), the MV-RPMD Hamiltonian is

HN=HRP+α=1N1βNxαTxα+1βNpαTpα1βNln|Θ|,
(5)

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

HRP=α=1NPαTPα2M+V0(Rα)+12MωN2(RαRα+1)T(RαRα+1),
(6)

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

Θ=Re(Tr[Γ]),
(7)

where

Γ=α=1N(Cα12I)M(Rα,Rα+1),
(8)
Cα=(xα+ipα)(xαipα)T,
(9)

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

Mnm(Rα,Rα+1)=eβN2[Vnn(Rα)+Vnn(Rα+1)]+O(βN2),n=mjnβN4[Vnj(Rα)+Vnj(Rα+1)]eβN2[Vjj(Rα)+Vjj(Rα+1)]  +jmβN4[Vjm(Rα)+Vjm(Rα+1)]eβN2[Vjj(Rα)+Vjj(Rα+1)]+O(βN2),nm,
(10)

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.

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 studies41 and surface hopping based simulations13 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,

H=Ps22ms+X,X,Y,Y=DA|XYVXYXY(s)XY|+jPj22M+12Mωj2(QjcjsMωj2)2,
(11)

where s, Ps, and ms are the position, momentum, and mass of the solvent polarization coordinate and VXY,XY(s) are the elements of the diabatic potential energy matrix where the subscripts X/Y/X/Y={D,A} label the donor and acceptor states of the particles. In Eq. (11), Pj, Qj, and M are the momentum, position, and mass of the jth bath mode, and cj is the coupling between the solvent and the jth bath mode of frequency ωj. The bath spectral density is Ohmic,

J(ω)=ηωeω/ωc,
(12)

with cut-off frequency ωc = ωs, and the dimensionless parameter η/msωs determines the coupling strength between the solvent and the bath modes.42 The continuous spectral density is discretized into f oscillators with frequencies23 

ωj=ωclogj0.5f,
(13)

and the coupling constants cj are defined as

cj=ωj2ηMωcfπ1/2,
(14)

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

VXYXY(s)=as2+bs+c,
(15)

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

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

CAB(t)=sgn(Θ)A({ξα}0)B({ξα}t)Wsgn(Θ)W,
(16)

where {ξα}t represents the set of bead positions and momenta {Rα, Pα, xα, pα} at time t, and the bead-averaged functions A({ξα}0)=1/NαA(ξα(0)) and B({ξα}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βNHN({ξα}0) with the MV-RPMD Hamiltonian, HN, 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,

Ṙα=HNPα,Ṗα=HNRα,ẋα=HNpα,ṗα=HNxα.
(17)

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,

CPn,h(t)=Trρneq(0)Pn(t)h,
(18)

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

h=h(sts),forwardh(sst),backward.
(19)

In the MV-RPMD framework, the Heaviside function in Eq. (18) is written in terms of the solvent ring polymer centroid, hh(±(s¯ts)), where s¯=1/Nα=1Nsα. The nth state populations at time t are evaluated using the “Boltzmann” estimator,31,34

Pnβ=ΓnnTr[Γ],
(20)

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, ρneq=ρneqsysρ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

ρneqsys=eβHsδ(ss)n=1Kδ(PnPn),
(21)

where Hs is the subsystem Hamiltonian given by the first line of Eq. (11), Pn is the population of the nth state, and the solvent position, s, and electron-proton state populations, Pn, 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

ρneqsys(0)=eβHRP(s)δ(s¯0s)n=1Kδ(PnSC(0)Pn),
(22)

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

PnSC=1Nα=1NPnSCα=12Nα=1N(xαn2+pαn21),
(23)

where PnSCα 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 = 039 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β to the same values at t = 0 since the latter includes the correct Boltzmann weights for each electronic state at a given nuclear configuration.

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.

TABLE I.

Solvent and bath parameters common to all three model PCET systems.

ParameteraModel IModel IIModel III
ms 22 000 22 000 22 000 
ωs × 104 3.72 4.00 3.72 
f 12 12 12 
M ms ms ms 
η/msωs 
T/K 300 300 300 
ParameteraModel IModel IIModel III
ms 22 000 22 000 22 000 
ωs × 104 3.72 4.00 3.72 
f 12 12 12 
M ms ms ms 
η/msωs 
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,

PnSCα=Pn.
(24)

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 = 0 a.u., and only the DD and AA states are populated with PDD=PAA=0.5 and PDA=PAD=0. For each model, we sample the distribution with a total of 5 × 108 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 × 104, 8 × 104, and 1.5 × 105 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 = −0.8 a.u. (at the intersection of the DD and AD states), and the initial electronic state populations are taken to be PDD=PAD=0.5 and PDA=PAA=0. For this simulation, trajectories were integrated for a total time of 1200 fs and 2.5 × 104 trajectories were employed to obtain the converged results shown here.

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.

FIG. 1.

The quasi-diabatic state potentials as a function of solvent coordinate are shown for model I, with state DD in red, DA in green, AD in blue, and AA in pink.

FIG. 1.

The quasi-diabatic state potentials as a function of solvent coordinate are shown for model I, with state DD in red, DA in green, AD in blue, and AA in pink.

Close modal
FIG. 2.

Population dynamics for model I (concerted), where population transfers directly from the reactant DD state (in red) to the product AA state in pink. The intermediate AD (in blue) and DA (in green) states are not populated during the course of the reaction.

FIG. 2.

Population dynamics for model I (concerted), where population transfers directly from the reactant DD state (in red) to the product AA state in pink. The intermediate AD (in blue) and DA (in green) states are not populated during the course of the reaction.

Close modal

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.

FIG. 3.

The quasi-diabatic state potentials as a function of solvent coordinate for model II with state DD in red, DA in green, AD in blue, and AA in pink.

FIG. 3.

The quasi-diabatic state potentials as a function of solvent coordinate for model II with state DD in red, DA in green, AD in blue, and AA in pink.

Close modal
FIG. 4.

Population dynamics for model II (sequential proton transfer followed by the electron transfer), where population first transfers from the reactant DD state (in red) to the DA state (in green) corresponding to the proton transfer before the electron transfers leading to a rapid rise in the population of the product AA state in pink. There is a small thermal population in the AD state in blue.

FIG. 4.

Population dynamics for model II (sequential proton transfer followed by the electron transfer), where population first transfers from the reactant DD state (in red) to the DA state (in green) corresponding to the proton transfer before the electron transfers leading to a rapid rise in the population of the product AA state in pink. There is a small thermal population in the AD state in blue.

Close modal

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.

FIG. 5.

The quasi-diabatic state potentials as a function of solvent coordinate are shown for model III, with state DD in red, DA in green, AD in blue, and AA in pink.

FIG. 5.

The quasi-diabatic state potentials as a function of solvent coordinate are shown for model III, with state DD in red, DA in green, AD in blue, and AA in pink.

Close modal
FIG. 6.

Population dynamics for model III (sequential ET-PT), where population first transfers from the reactant DD state (in red) to the AD state (in blue) corresponding to the electron transfer before the proton transfers leading to a rapid rise in the population of the product AA state in pink. The DA state (in green) shows some initial thermal population but is not populated during the course of the reaction.

FIG. 6.

Population dynamics for model III (sequential ET-PT), where population first transfers from the reactant DD state (in red) to the AD state (in blue) corresponding to the electron transfer before the proton transfers leading to a rapid rise in the population of the product AA state in pink. The DA state (in green) shows some initial thermal population but is not populated during the course of the reaction.

Close modal

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.

FIG. 7.

Population dynamics for model III (sequential ET-PT), with the reactant state in red, PT state in green, ET state in blue, and product state in pink where trajectories are initialized to the electron transfer transition state.

FIG. 7.

Population dynamics for model III (sequential ET-PT), with the reactant state in red, PT state in green, ET state in blue, and product state in pink where trajectories are initialized to the electron transfer transition state.

Close modal

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 theory45 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

kFGR=2πω|Δ|2evzScothzIvS cschz,
(25)

where ω=2a/ms is the frequency of the product diabatic state, v = (VRVP)/ω, z = βω/2, S=msωVd2/2, Iv is a modified Bessel function of the first kind, Vd is the horizontal displacement of the diabatic potential energy functions, and VR/P are the values of the potential energy at the reactant/product minimum such that VRVP measures the driving force. For adiabatic ET, we use Kramers theory,45 

kKT=1+γ2ωb2γ2ωbω2πeβGcl,
(26)

where ωb is the frequency at the top of the barrier, Gcl is the solvent FE barrier when the solvent is treated classically, and γ = η/MS.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.

TABLE II.

FGR and Kramers theory rates (indicated with a *) for the concerted PCET (kDD→AA), electron transfer (kDD→AD), and proton transfer (kDD→DA) from the reactant DD state for all three models are reported in s−1. The fastest rate for each model is highlighted in bold to indicate the preferred mechanism.

Reaction pathModel IModel IIModel III
kDD→AA 1.85 × 107 1.61 × 106 4.70 × 106 
kDD→DA 9.81 × 10−17 2.53 × 109 5.97 × 104 
kDD→AD 2.69 × 105* 1.01 × 106 1.03 × 1011* 
Reaction pathModel IModel IIModel III
kDD→AA 1.85 × 107 1.61 × 106 4.70 × 106 
kDD→DA 9.81 × 10−17 2.53 × 109 5.97 × 104 
kDD→AD 2.69 × 105* 1.01 × 106 1.03 × 1011* 

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.

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.

In the limit that N, the high-temperature symmetric Trotter approximation is used to separate the state independent nuclear potential operator, V0, and the diabatic potential energy matrix, V, from the nuclear kinetic energy operator T,

n,Rα|eβNH|Rα+1,mn,Rα|eβN2V0eβN2VeβNTeβN2VeβN2V0|Rα+1,m=eβN2(V0(Rα)+V0(Rα+1))Rα|eβNT|Rα+1  ×n|eβN2V(Rα)eβN2V(Rα+1)|m.
(A1)

The nuclear kinetic matrix element can be evaluated exactly to obtain

Rα|eβNT|Rα+1=dPRα|PP|eβNT|Rα+1=dPRα|PeβNP2/2mP|Rα+1=M2πβN1/2eβN2MωN2(RαRα+1)2.
(A2)

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

n,Rα|eβNHN|Rα+1,mM2πβN1/2eβN2V0Rα+V0(Rα+1)+MωN22(RαRα+1)2  ×n|eβN2[V(Rα)+V(Rα+1)]|m.
(A3)

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

n|eβN(VD+VOD)|mn|eβN2VODeβNVDeβN2VOD|m=j,kn|eβN2VOD|jj|eβNVD|kk|eβN2VOD|m=jn|eβN2VOD|jeβNVjjj|eβN2VOD|m.
(A4)

The off-diagonal matrix elements are easily evaluated,

n|eβN2VOD|jn|(1βN2VOD)|j+O(βN2)=1,n=jβN2Vnj,nj
(A5)

where Vnj 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,

n|eβN(VD+VOD)|n=n|eβN2VOD|neβNVnnn|eβN2VOD|n=eβNVnn.
(A6)
  1. If nj,

n|eβN(VD+VOD)|n=njn|eβN2VOD|jeβNVjjj|eβN2VOD|n=njβN22Vnj2eβNVjj=0+O(βN2).
(A7)
  • (nm):

    If n = j and if mj,

mjn|eβN2VOD|jeβNVjjj|eβN2VOD|m=mjβN2VjmeβNVjj.
(A8)
  1. If nj and if m = j,

njn|eβN2VOD|jeβNVjjj|eβN2VOD|m=njβN2VnjeβNVjj.
(A9)
  1. If nj and if mj,

nmjn|eβN2VOD|jeβNVjjj|eβN2VOD|m=nmjβN22VnjVjmeβNVjj=0+O(βN2).
(A10)

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

H=Ps22ms+PR22mR+Vp(R)+Vps(R,s)+Vij(R,s).
(B1)

In Eq. (B1), R is the proton coordinate with conjugate momentum PR, and Vp(R) is a double well potential in the proton coordinate,

Vp(R)=mRωR22R2+mR2ωR416V0R4λR3,
(B2)

where mR is the mass of the proton, ωR is the frequency, λ is a measure of anharmonicity, and V0 determines the height of the barrier for the proton transfer. Further, the proton-solvent coupling is

Vps(R,s)=μ1stanh(ϕR),
(B3)

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

Vii(R,s)=12msωs2(ssi)2+aiμ2tanh(ϕR),
(B4)

where μ2, ai, 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.

TABLE III.

Parameters for the model Hamiltonians in Eq. (B1).

ParameteraModel IModel IIModel III
mR 1836.1 1836.1 1836.1 
ωR 0.0104 0.0104 0.0104 
V0 0.012 0.014 0.012 
s1 −2.13 −2.16 −2.13 
s2 2.13 2.16 2.13 
V12 0.002 45 0.0124 0.002 45 
μ1 0.001 1 0.017 −0.001 1 
μ2 × 103 5.84 0.71 5.84 
λ 0.0 0.012 0.0 
ParameteraModel IModel IIModel III
mR 1836.1 1836.1 1836.1 
ωR 0.0104 0.0104 0.0104 
V0 0.012 0.014 0.012 
s1 −2.13 −2.16 −2.13 
s2 2.13 2.16 2.13 
V12 0.002 45 0.0124 0.002 45 
μ1 0.001 1 0.017 −0.001 1 
μ2 × 103 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 −6a0s ≤ 6a0, we diagonalize the system Hamiltonian on a uniform Discrete Variable Representation (DVR) grid in the proton coordinate with a grid range of −2a0R ≤ 2a0 and 100 grid points. The adiabatic eigenstates obtained upon diagonalizing the system Hamiltonian are written as R;s|ϵi, where ϵi is the ith adiabatic state with eigenenergy Ei.

Further, by diagonalizing the system Hamiltonian for a single electronic state (donor or acceptor) at each value of s, we construct localized proton wavefunctions, R;s|lj, where lj is the jth quasi-diabatic local electron-proton states that can be expressed in terms of the adiabatic eigenstates as

R;s|lj=idRR;s|ϵiϵi|R;sR;s|lj.
(B5)

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

lj|H|lj=i,ilj|ϵiϵi|H|ϵiϵi|lj=ilj|ϵiEiϵi|lj,
(B6)

where Ei is the energy of the ith 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

ϵi|lj=dRϵi|RR|lj.
(B7)

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

TABLE IV.

Diabatic potential energy surface parameters for model I.

Diabaticabc
VDD 0.0015 0.0075 −0.0041 
VDA 0.0015 0.0055 0.0072 
VAD 0.0015 −0.0055 0.0072 
VAA 0.0015 −0.0075 −0.0041 
Diabaticabc
VDD 0.0015 0.0075 −0.0041 
VDA 0.0015 0.0055 0.0072 
VAD 0.0015 −0.0055 0.0072 
VAA 0.0015 −0.0075 −0.0041 
TABLE V.

Diabatic coupling matrix elements for model I.

CouplingΔ
VDD,DA 9.7 × 10−5 
VDD,AD 2.5 × 10−3 
VDD,AA 1.8 × 10−4 
VDA,AD 1.8 × 10−4 
VDA,AA 2.5 × 10−3 
VAD,AA 9.7 × 10−5 
CouplingΔ
VDD,DA 9.7 × 10−5 
VDD,AD 2.5 × 10−3 
VDD,AA 1.8 × 10−4 
VDA,AD 1.8 × 10−4 
VDA,AA 2.5 × 10−3 
VAD,AA 9.7 × 10−5 
TABLE VI.

Diabatic potential energy surface parameters for model II.

Diabaticabc
VDD 0.0015 0.0072 −0.0018 
VDA 0.0018 0.0058 −0.0013 
VAD 0.0018 −0.0061 0.0034 
VAA 0.0016 −0.0083 −0.0018 
Diabaticabc
VDD 0.0015 0.0072 −0.0018 
VDA 0.0018 0.0058 −0.0013 
VAD 0.0018 −0.0061 0.0034 
VAA 0.0016 −0.0083 −0.0018 
TABLE VII.

Diabatic coupling matrix elements for model II.

CouplingΔ
VDD,DA 1.1 × 10−3 
VDD,AD 1.2 × 10−4 
VDD,AA 1.2 × 10−4 
VDA,AD 1.2 × 10−4 
VDA,AA 1.2 × 10−4 
VAD,AA 1.4 × 10−3 
CouplingΔ
VDD,DA 1.1 × 10−3 
VDD,AD 1.2 × 10−4 
VDD,AA 1.2 × 10−4 
VDA,AD 1.2 × 10−4 
VDA,AA 1.2 × 10−4 
VAD,AA 1.4 × 10−3 
TABLE VIII.

Diabatic potential energy surface parameters for model III.

Diabaticabc
VDD 0.0015 0.008 0.0009 
VDA 0.0015 0.0098 0.013 
VAD 0.0015 −0.0056 −0.0095 
VAA 0.0015 −0.013 0.0009 
Diabaticabc
VDD 0.0015 0.008 0.0009 
VDA 0.0015 0.0098 0.013 
VAD 0.0015 −0.0056 −0.0095 
VAA 0.0015 −0.013 0.0009 
TABLE IX.

Diabatic coupling matrix elements for model III.

CouplingΔ
VDD,DA 6.9 × 10−4 
VDD,AD 2.5 × 10−3 
VDD,AA 1.8 × 10−4 
VDA,AD 1.8 × 10−4 
VDA,AA 2.5 × 10−3 
VAD,AA 6.9 × 10−4 
CouplingΔ
VDD,DA 6.9 × 10−4 
VDD,AD 2.5 × 10−3 
VDD,AA 1.8 × 10−4 
VDA,AD 1.8 × 10−4 
VDA,AA 2.5 × 10−3 
VAD,AA 6.9 × 10−4 
1.
M. H. V.
Huynh
and
T. J.
Meyer
,
Chem. Rev.
107
,
5004
(
2007
).
2.
R. I.
Cukier
and
D. G.
Nocera
,
Annu. Rev. Phys. Chem.
49
,
337
(
1998
).
3.
S.
Hammes-Schiffer
and
A. V.
Soudackov
,
J. Phys. Chem. B
112
,
14108
(
2008
).
4.
J. J.
Warren
,
T. A.
Tronic
, and
J. M.
Mayer
,
Chem. Rev.
110
,
6961
(
2010
).
7.
S.
Hammes-Schiffer
and
A. A.
Stuchebrukhov
,
Chem. Rev.
110
,
6939
(
2010
).
8.
S.
Hammes-Schiffer
,
J. Am. Chem. Soc.
137
,
8860
(
2015
).
9.
A. V.
Soudackov
and
S.
Hammes-Schiffer
,
Faraday Discuss.
195
,
171
(
2016
).
10.
M. T.
Huynh
,
S. J.
Mora
,
M.
Villalba
,
M. E.
Tejeda-Ferrari
,
P. A.
Liddell
,
B. R.
Cherry
,
A.-L.
Teillout
,
C. W.
Machan
,
C. P.
Kubiak
,
D.
Gust
,
T. A.
Moore
,
S.
Hammes-Schiffer
, and
A. L.
Moore
,
ACS Cent. Sci.
3
,
372
(
2017
).
11.
P.
Goyal
,
C. A.
Schwerdtfeger
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. B
119
,
2758
(
2015
).
12.
P.
Goyal
,
C. A.
Schwerdtfeger
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. B
120
,
2407
(
2016
).
13.
J.-Y.
Fang
and
S.
Hammes-Schiffer
,
J. Chem. Phys.
106
,
8442
(
1997
).
14.
H.
Decornez
and
S.
Hammes-Schiffer
,
J. Phys. Chem. A
104
,
9370
(
2000
).
15.
A.
Soudackov
,
E.
Hatcher
, and
S.
Hammes-Schiffer
,
J. Chem. Phys.
122
,
014505
(
2004
).
16.
E.
Hatcher
,
A.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. B
109
,
18565
(
2005
).
17.
A.
Hazra
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. B
114
,
12319
(
2010
).
18.
A.
Hazra
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. Lett.
2
,
36
(
2011
).
19.
A.
Warshel
,
J. Phys. Chem.
86
,
2218
(
1982
).
20.
S. J.
Cotton
,
K.
Igumenshchev
, and
W. H.
Miller
,
J. Chem. Phys.
141
,
084104
(
2014
).
21.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
1965
).
22.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
 III
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
23.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
123
,
034102
(
2005
).
24.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
122
,
084106
(
2005
).
25.
T. J. H.
Hele
and
S. C.
Althorpe
,
J. Chem. Phys.
138
,
084108
(
2013
).
26.
S. C.
Althorpe
and
T. J. H.
Hele
,
J. Chem. Phys.
139
,
084115
(
2013
).
27.
J. S.
Kretchmer
and
T. F.
Miller
 III
,
J. Chem. Phys.
138
,
134109
(
2013
).
28.
J. S.
Kretchmer
and
T. F.
Miller
 III
,
Inorg. Chem.
55
,
1022
(
2016
).
29.
A. R.
Menzeleev
,
N.
Ananth
, and
T. F.
Miller
 III
,
J. Chem. Phys.
135
,
074106
(
2011
).
30.
P.
Shushkov
,
R.
Li
, and
J. C.
Tully
,
J. Chem. Phys.
137
,
22A549
(
2012
).
31.
N.
Ananth
,
J. Chem. Phys.
139
,
124102
(
2013
).
32.
J. O.
Richardson
and
M.
Thoss
,
J. Chem. Phys.
139
,
031102
(
2013
).
33.
A. R.
Menzeleev
,
F.
Bell
, and
T. F.
Miller
 III
,
J. Chem. Phys.
140
,
064103
(
2014
).
34.
J. R.
Duke
and
N.
Ananth
,
J. Phys. Chem. Lett.
6
,
4219
(
2015
).
35.
H.-D.
Meyer
and
W. H.
Miller
,
J. Chem. Phys.
70
,
3214
(
1979
).
36.
G.
Stock
and
M.
Thoss
,
Phys. Rev. Lett.
78
,
578
(
1997
).
37.
N.
Ananth
and
T. F.
Miller
 III
,
J. Chem. Phys.
133
,
234103
(
2010
).
38.
H. F.
Trotter
,
Proc. Am. Math. Soc.
10
,
545
(
1958
).
39.
T. J. H.
Hele
and
N.
Ananth
,
Faraday Discuss.
195
,
269
(
2016
).
40.
D.
Chandler
,
Introduction to Modern Statistical Mechanics
(
Oxford University Press
,
1987
).
41.
N.
Ananth
and
T. F.
Miller
 III
,
Mol. Phys.
110
,
1009
(
2012
).
42.
A.
Caldeira
and
A.
Leggett
,
Ann. Phys.
149
,
374
(
1983
).
43.
J. O.
Richardson
,
P.
Meyer
,
M.-O.
Pleinert
, and
M.
Thoss
,
Chem. Phys.
482
,
124
(
2017
).
44.
E.
Fermi
,
Nuclear Physics: A Course Given by Enrico Fermi at the University of Chicago
(
University of Chicago Press
,
1950
).
45.
N. E.
Henriksen
and
F. Y.
Hansen
,
Theories of Molecular Reaction Dynamics: The Microscopic Foundation of Chemical Kinetics
(
Oxford University Press on Demand
,
2008
).
46.
J.
Ulstrup
and
J.
Jortner
,
J. Chem. Phys.
63
,
4358
(
1975
).
47.
J.
Ulstrup
,
Charge Transfer Processes in Condensed Media
(
Springer Science & Business Media
,
2012
).
48.
N.
Makri
,
J. Chem. Phys.
106
,
2286
(
1997
).
49.
R. I.
Cukier
,
J. Phys. Chem.
98
,
2377
(
1994
).
50.
R. I.
Cukier
,
J. Phys. Chem.
99
,
16101
(
1995
).
51.
R. I.
Cukier
,
J. Phys. Chem.
100
,
15428
(
1996
).