Rate constant expressions for vibronically nonadiabatic proton transfer and proton-coupled electron transfer reactions are presented and analyzed. The regimes covered include electronically adiabatic and nonadiabatic reactions, as well as high-frequency and low-frequency proton donor-acceptor vibrational modes. These rate constants differ from previous rate constants derived with the cumulant expansion approach in that the logarithmic expansion of the vibronic coupling in terms of the proton donor-acceptor distance includes a quadratic as well as a linear term. The analysis illustrates that inclusion of this quadratic term in the framework of the cumulant expansion framework may significantly impact the rate constants at high temperatures for proton transfer interfaces with soft proton donor-acceptor modes that are associated with small force constants and weak hydrogen bonds. The effects of the quadratic term may also become significant in these regimes when using the vibronic coupling expansion in conjunction with a thermal averaging procedure for calculating the rate constant. In this case, however, the expansion of the coupling can be avoided entirely by calculating the couplings explicitly for the range of proton donor-acceptor distances sampled. The effects of the quadratic term for weak hydrogen-bonding systems are less significant for more physically realistic models that prevent the sampling of unphysical short proton donor-acceptor distances. Additionally, the rigorous relation between the cumulant expansion and thermal averaging approaches is clarified. In particular, the cumulant expansion rate constant includes effects from dynamical interference between the proton donor-acceptor and solvent motions and becomes equivalent to the thermally averaged rate constant when these dynamical effects are neglected. This analysis identifies the regimes in which each rate constant expression is valid and thus will be important for future applications to proton transfer and proton-coupled electron transfer in chemical and biological processes.

Vibronic couplings play a key role in determining the rate constants for proton transfer (PT) and proton-coupled electron transfer (PCET). In the Golden Rule formulation, the rate constant is proportional to the square of the vibronic coupling at the environmental configurations corresponding to the degeneracy of the reactant and product vibronic states. In this context, the vibronic coupling is the Hamiltonian matrix element between the reactant and product diabatic electron-proton vibronic states. Typically, these vibronic couplings are assumed to be only weakly dependent on the configuration of the environment (i.e., the Condon approximation). In the case of PT or PCET reactions, however, the proton donor-acceptor distance coordinate represents an exception. In typical two-dimensional (2D) potentials describing the PT interface in terms of the proton coordinate q and the proton donor-acceptor distance R, the coordinates q and R are strongly coupled, resulting in the strong dependence of the PT barrier height along the proton coordinate q on the proton donor-acceptor distance R. A major consequence of this effect is the strong dependence of the vibronic coupling on the proton donor-acceptor distance.

The importance of the dependence of the vibronic coupling on the proton donor-acceptor vibrational motion has been recognized in the literature. The effects of the proton donor-acceptor motion were explicitly included in theoretical treaments of PT1–5 and PCET6–8 reactions. The dependence of the nonadiabatic vibronic coupling V on the proton donor-acceptor distance R is often represented in exponential form by expanding the natural logarithm of the coupling around the equilibrium proton donor-acceptor distance R0. When only the terms up to second order in δR = RR0 are retained, this dependence is expressed as

lnV(R)=lnV(R0)αδR12γδR2,
(1)

where

α=dlnV(R)dRR=R0=1V(R0)dV(R)dRR=R0
(2)

and

γ=d2lnV(R)dR2R=R0=1V(R0)d2V(R)dR2R=R0+1[V(R0)]2dV(R)dRR=R02.
(3)

Using this expansion, the vibronic coupling can be written in the following exponential form with V0V(R0):

V(R)=V0expαδR12γδR2.
(4)

Several previous treatments1,2,4,5,8,9 assumed that the terms linear in δR are dominant and that the quadratic term can be neglected. This approximation was justified by agreement with the full quantum mechanical treatment of both the q and R coordinates in various model PCET systems for the electronically nonadiabatic case, where the vibronic coupling is proportional to the R-dependent overlap integral between the reactant and product proton vibrational wavefunctions.8,10 In other cases,3,11,12 however, the quadratic term in Eq. (4) could potentially have a significant effect on the rate constant, especially in the case of the so-called deep tunneling or electronically adiabatic but vibronically nonadiabatic regime, where the vibronic coupling is represented by half the tunneling splitting in a double-well PT potential. These effects are also relevant to thermal averaging treatments of PT and PCET, in which the rate constant is averaged over the equilibrium distribution of proton donor-acceptor distances.6 In this treatment, the classical distribution function over the proton donor-acceptor distances may enable sampling of the relatively short distances at which the linear approximation for the logarithm of the vibronic coupling is insufficient, and the quadratic term may become necessary.

In this paper, we systematically analyze the effects of the quadratic term on the nonadiabatic rate constants for PT and PCET in various regimes. The aim of this study is to obtain deeper insight into how the characteristics of the PT interface, including the features of the potential energy surface and the effective mass and frequency of the proton donor-acceptor vibrational mode, impact the significance of the quadratic term. Section II summarizes the theoretical treatments for nonadiabatic PT and PCET reactions and provides expressions for the vibronic coupling in the relevant regimes. Section III presents results for a range of model systems, and Section IV summarizes the conclusions of this study.

The nonadiabatic theories of PT and PCET are based on the general Golden Rule formulation. In the conventional formulation, nonadiabatic transitions occur between the reactant and product quantum states interacting with the fluctuating environment in thermal equilibrium at a given temperature. The thermal fluctuations of the environment bring a pair of quantum states into degeneracy, thereby enabling a quantum transition with a probability that is proportional to the square of the quantum mechanical coupling between the corresponding states. Defining the reactant (I) and product (II) channels with Hamiltonians HˆI and HˆII, which give rise to the manifolds of reactant and product vibronic states for the entire system, the Golden Rule rate constant expression can be written as a thermodynamic average over the reactant states (i.e., eigenstates of the Hamiltonian HˆI) of the entire system:8 

k=1ħ2dtVˆ(0)exp()iħ0tΔHˆ(τ)dτVˆ(t){I}=1ħ2J(t)dt.
(5)

The time-dependent quantities in the above expression are op erators for the coupling Vˆ and the energy gap ΔHˆHˆIIHˆI in the Heisenberg representation:

(6)
Vˆ(t)=eiħHˆItVˆeiħHˆIt,
ΔHˆ(t)=eiħHˆItHˆIIHˆIeiħHˆIt.
The exp(−) denotes a negatively time-ordered exponential, and {I}TrexpβHˆI/TrexpβHˆI, where β−1 = kBT, denotes a quantum mechanical thermal average over initial states.

The reactant and product vibronic states form two sets of quantum states in the electronically diabatic representation. As in Marcus-Levich-Dogonadze nonadiabatic ET theory,13,14 the diabatic electronic states are characterized by the reactant and product electronic charge distributions in the reaction complex. These diabatic electronic states are defined either by the transferring electron localized on the electron donor or acceptor species, as in the case of ET or PCET, or by the specific bonding pattern describing the proton bonded to the proton donor or acceptor, as in the case of PT. For PCET and PT, the proton donor-acceptor distance R warrants special treatment due to its special role in modulating the couplings between the reactant and product vibronic states. In the remainder of this section, we consider different treatments of this coordinate, emphasizing the effects of the linear and quadratic terms in the logarithmic expansion of the coupling matrix element given in Eq. (4). In the general rate constant expressions, we do not specify the form of the vibronic coupling, assuming that any coupling matrix element can be represented in the form of Eq. (4) near the equilibrium value of the proton donor-acceptor distance. Specific expressions for the vibronic coupling in the electronically adiabatic and nonadiabatic limits, as well as the intermediate regime, are provided at the end of Section II A.

The reactant and product vibronic states are the eigenstates of the reactant and product channel Hamiltonians HˆI and HˆII, respectively. Within each of these channels, the reactant and product states can be represented in the following form:

ΨμmJ(re,q,R,ξ)=ΦμJ(re,q|R,ξ)χμmJ(R,ξ),J=I,II.
(7)

Here, ΦμJ are the diabatic electron-proton vibronic wavefunctions parametrically dependent on the donor-acceptor distance R and the environmental coordinates ξ for the reactant (J = I) and product (J = II) vibronic manifolds with corresponding diabatic electronic charge distributions. The corresponding vibronic energies, EμJ(R,ξ), serve as vibronic potential energy surfaces for the remaining degrees of freedom, and the final quantization along R and ξ provides the vibrational wavefunctions of the environment, χμmJ(R,ξ). In the linear response approximation, the vibronic surfaces EμJ(R,ξ) are represented as harmonic surfaces in the multidimensional space spanned by the proton donor-acceptor coordinate R and the environmental modes.

The matrix elements of the coupling operator Vˆ between the reactant and product states defined in Eq. (7) are given by

ΨμmIVˆΨνnII=χμmIVμν(R,ξ)χνnII,
(8)

where

Vμν(R,ξ)=ΦμIVˆΦνIIre,qVμν(R)
(9)

is the electron-proton vibronic coupling appearing in Eq. (4). According to the Condon approximation, this coupling is assumed to be independent of the environmental coordinates, but its dependence on the proton donor-acceptor distance cannot be neglected.

In the electronically diabatic channel representation, the hierarchy of adiabatic approximations between the electrons and the proton, as well as between the proton and the slower nuclear degrees of freedom within each channel, is usually well-justified. The adiabatic separation between the proton coordinate q and the proton donor-acceptor coordinate R is less obvious because these coordinates are strongly coupled in typical two-dimensional potentials describing hydrogen-bonded PT interfaces. However, previous work15 has established that this adiabatic separation, in which the proton coordinate q is treated as a “fast” coordinate and the proton donor-acceptor coordinate R is treated as a “slow” coordinate, is well-justified not only for electronically diabatic but also for a wide range of electronically adiabatic potentials, including double well potentials with high barriers (i.e., in the deep tunneling regime).

The vibronic coupling V(R) can be calculated from the proton potential energy curves associated with the reactant and product diabatic electronic states, as depicted in Figure 1(a), in conjunction with the electronic coupling Vel between these two states. Georgievskii and Stuchebrukhov16 derived a semiclassical expression for the general vibronic coupling V(sc),

V(sc)=κV(ad),
(10)

where V(ad) = Δ(tun)/2 is the vibronic coupling in the electronically adiabatic limit, equal to the half the tunneling splitting between the proton vibrational levels in the ground state adiabatic electronic potential (Figure 1(b)), and

κ=2πp12eplnppΓp+1.
(11)

In this expression, Γ(x) is the gamma function, and p = τp/τe is the electron-proton adiabaticity parameter defined in terms of the effective electronic transition time,

τe=ħVel,
(12)

and the effective proton tunneling time,

τp=VelΔFvt.
(13)

Here, Vel is the electronic coupling between the reactant and product diabatic electronic states,17F| is the difference between the slopes of the electronically diabatic potentials at the crossing point, and the tunneling velocity is

vt=2VEm,
(14)

where V is the energy at which the potential energy curves cross, E is the tunneling energy associated with the unperturbed degenerate proton vibrational levels in the reactant and product diabatic potentials, and m is the mass of the proton.

FIG. 1.

(a) Reactant (blue) and product (red) electronically diabatic proton potentials and associated proton vibrational wavefunctions at a fixed proton donor acceptor distance; (b) ground state electronically adiabatic proton potential and the corresponding ground (blue) and excited (red) state proton vibrational wavefunctions at a fixed proton donor-acceptor distance. The tunneling splitting is the energy difference between these two proton vibrational energy levels, as indicated by dashed lines.

FIG. 1.

(a) Reactant (blue) and product (red) electronically diabatic proton potentials and associated proton vibrational wavefunctions at a fixed proton donor acceptor distance; (b) ground state electronically adiabatic proton potential and the corresponding ground (blue) and excited (red) state proton vibrational wavefunctions at a fixed proton donor-acceptor distance. The tunneling splitting is the energy difference between these two proton vibrational energy levels, as indicated by dashed lines.

Close modal

The semiclassical expression for the vibronic coupling given in Eq. (10) spans the electronically adiabatic and nonadiabatic regimes. In the electronically adiabatic limit, p ≫ 1, κ = 1, and the vibronic coupling simplifies to V(ad). Several semiclassical expressions in terms of parameters associated with proton double well potentials have been derived for V(ad).5,11,12,18,19 In the electronically nonadiabatic limit, p ≪ 1, κ=2πp1/2, and the vibronic coupling reduces to

V(na)=Velφ(I)φ(II),
(15)

where φ(I) and φ(II) are the proton vibrational wavefunctions in the reactant and product wells, respectively (Figure 1(a)).

From a physical perspective, the electronically adia batic limit corresponds to the case in which the electrons respond instantaneously to the proton motion, resulting in the reaction proceeding entirely on the ground adiabatic electronic state. In contrast, the electronically nonadiabatic limit corresponds to the case in which the electrons are unable to respond instantaneously to the proton motion, thereby leading to involvement of the excited electronic state. Most PT and hydrogen atom transfer (HAT) reactions are electronically adiabatic, occurring on the ground adiabatic electronic state. In contrast, PCET reactions that involve concerted electron and proton transfer between different donors and acceptors (i.e., different orbitals, atoms, or groups) are typically electronically nonadiabatic. The distinctions between these types of reactions have been discussed extensively elsewhere.20–24 

In the approaches used previously for PT2,5 and PCET8 reactions in solution, the proton donor-acceptor coordinate R is assumed to be adiabatically separated from the electron and proton coordinates. Here, we limit our discussion to a single pair of reactant and product electron-proton vibronic states ΦμI and ΦνII to avoid complicated notation; the total rate constant is obtained by averaging over the reactant states and summing over the product states.

Expanding the energy gap in a Taylor series in the coordinate R up to the linear term and substituting the expansion in Eq. (4) for the vibronic coupling, we obtain the following expression for the reaction flux correlation function, J(t), in Eq. (5):5 

J(t)=V02expiħEˆteiħα̃δRˆ(0)i2ħγ̃δRˆ2(0)exp()iħ0tδEˆ(τ)+Λ̃RδRˆ(τ)dτeiħα̃δRˆ(t)i2ħγ̃δRˆ2(t){I},
(16)

where δEˆ(τ) and δRˆ(τ) are the Heisenberg operators describing the time evolution of fluctuations of the energy gap EˆΔHˆ and proton donor-acceptor distance R, respectively. Moreover, α̃=iħα, γ̃=iħγ, and Λ̃R is the derivative of the energy gap with respect to the proton donor-acceptor coordinate R evaluated at the minimum of the reactant vibronic energy surface. In the harmonic approximation, Λ̃R is proportional to ΔR, the shift of the equilibrium value of the proton donor-acceptor distance in the product state relative to that in the reactant state. Note that for symmetric systems, the reactant and product equilibrium proton donor-acceptor distances are equal, and thus, Λ̃R=0.

Expanding the exponentials in Eq. (16) and collecting terms of the same order in i/ℏ, we obtain a series that can be represented as a cumulant expansion of the time-ordered exponential for the time-dependent variable Xˆ(τ),

exp()iħ0tXˆ(τ)dτ=expiħ0tdτXˆ(τ)+iħ20tdτ10τ1dτ2Xˆ(τ2)Xˆ(τ1)+,
(17)

where

Xˆ(τ)=α̃tδRˆ(0)+γ̃2tδRˆ2(0)+δEˆ(τ)+Λ̃RδRˆ(τ)+α̃tδRˆ(t)+γ̃2tδRˆ2(t),0<τ<t
(18)

and double angular brackets denote the cumulants of Xˆ(τ), which can be expressed in terms of the central moments and correlation functions of the time-dependent fluctuations of the energy gap and proton donor-acceptor distance, δEˆ and δRˆ, respectively. These moments can in turn be expressed in terms of the corresponding cumulants. Next, we assume that these time-dependent fluctuations can be described by Gaussian processes corresponding to the time evolution on harmonic surfaces. For the proton donor-acceptor coordinate R, the Gaussian process δRˆ(t) corresponds to the time evolution of a harmonic oscillator with an effective mass M and frequency Ω. In practice, the effective mass and frequency can be defined in terms of the normal modes of the system that have nonvanishing projections on the proton donor-acceptor axis and thus contribute to the fluctuations of the coordinate R.25 In addition, we assume that the corresponding harmonic potential along R is uncoupled from the other environmental degrees of freedom, which can be achieved by considering the normal modes of the entire system in the reactant state.

For Gaussian processes δEˆ(τ) and δRˆ(τ), all of the corresponding cumulants higher than second order vanish exactly. However, the cumulant expansion for the rate constant based on the expansion of the time-ordered exponential with the variable Xˆ(τ) does not truncate even for Gaussian δEˆ(τ) and δRˆ(τ) due to the presence of the quadratic terms δRˆ2(t) in Eq. (18). The final expression for the reaction flux can be written as

J(t)=V02expiħEtexpα2δR2+CR(t)2iαħΛ̃R0tCR(τ)dτ+n=11nγn1nKn(1)(t)+α2Kn+1(2)(t)1ħ20tdτ10τ1dτ2CE(τ1τ2)Λ̃R2ħ20tdτ10τ1dτ2CR(τ1τ2),
(19)

where

E=ΔG+Λ
(20)

and the operator notation is omitted for simplicity. Here, ΔG and Λ are the reaction free energy and the total reorganization energy, respectively, and the correlation functions are defined as

(21)
CE(τ1τ2)=δE(τ2)δE(τ1),
CR(τ1τ2)=δR(τ2)δR(τ1).
Eq. (19) is similar to the analogous expression derived in Ref. 5 for PT and Ref. 8 for PCET reactions except for an additional time-dependent term with Kn(1)(t) and Kn(2)(t) denoting nth order polynomials in CR(t) and δR2=CR(0). This additional term arises because of the quadratic term in the logarithmic expansion of the coupling. The explicit expressions for Kn(1)(t) and Kn(2)(t) for n = 1–8 (i.e., up to eighth order in γ) are given in Table I.

TABLE I.

Explicit expressions for polynomials Kn(1)(t) and Kn+1(2)(t) in Eq. (19) for n = 1–8 (up to eighth order in γ). The polynomials are given in terms of cumulants κ2,0δR2 and κ1,1δR(0)δR(t).

nKn(1)(t)Kn+1(2)(t)
κ2,0 κ1,12+2κ2,0κ1,1+κ2,02 
κ1,12+κ2,02 κ1,13+3κ2,0κ1,12+3κ2,02κ1,1+κ2,03 
κ2,03+3κ1,12κ2,0 κ1,14+4κ2,0κ1,13+6κ2,02κ1,12+4κ2,03κ1,1+κ2,04 
κ1,14+6κ2,02κ1,12+κ2,04 κ1,15+5κ2,0κ1,14+10κ2,02κ1,13+10κ2,03κ1,12+5κ2,04κ1,1+κ2,05 
κ2,05+10κ1,12κ2,03+5κ1,14κ2,0 κ1,16+6κ2,0κ1,15+15κ2,02κ1,14+20κ2,03κ1,13+15κ2,04κ1,12+6κ2,05κ1,1+κ2,06 
κ1,16+15κ2,02κ1,14+15κ2,04κ1,12+κ2,06 κ1,17+7κ2,0κ1,16+21κ2,02κ1,15+35κ2,03κ1,14+35κ2,04κ1,13+21κ2,05κ1,12+7κ2,06κ1,1+κ2,07 
κ2,07+21κ1,12κ2,05+35κ1,14κ2,03+7κ1,16κ2,0 κ1,18+8κ2,0κ1,17+28κ2,02κ1,16+56κ2,03κ1,15+70κ2,04κ1,14+56κ2,05κ1,13+28κ2,06κ1,12+8κ2,07κ1,1+κ2,08 
κ1,18+28κ2,02κ1,16+70κ2,04κ1,14+28κ2,06κ1,12+κ2,08 κ1,19+9κ2,0κ1,18+36κ2,02κ1,17+84κ2,03κ1,16+126κ2,04κ1,15+126κ2,05κ1,14+84κ2,06κ1,13+36κ2,07κ1,12+9κ2,08κ1,1+κ2,09 
nKn(1)(t)Kn+1(2)(t)
κ2,0 κ1,12+2κ2,0κ1,1+κ2,02 
κ1,12+κ2,02 κ1,13+3κ2,0κ1,12+3κ2,02κ1,1+κ2,03 
κ2,03+3κ1,12κ2,0 κ1,14+4κ2,0κ1,13+6κ2,02κ1,12+4κ2,03κ1,1+κ2,04 
κ1,14+6κ2,02κ1,12+κ2,04 κ1,15+5κ2,0κ1,14+10κ2,02κ1,13+10κ2,03κ1,12+5κ2,04κ1,1+κ2,05 
κ2,05+10κ1,12κ2,03+5κ1,14κ2,0 κ1,16+6κ2,0κ1,15+15κ2,02κ1,14+20κ2,03κ1,13+15κ2,04κ1,12+6κ2,05κ1,1+κ2,06 
κ1,16+15κ2,02κ1,14+15κ2,04κ1,12+κ2,06 κ1,17+7κ2,0κ1,16+21κ2,02κ1,15+35κ2,03κ1,14+35κ2,04κ1,13+21κ2,05κ1,12+7κ2,06κ1,1+κ2,07 
κ2,07+21κ1,12κ2,05+35κ1,14κ2,03+7κ1,16κ2,0 κ1,18+8κ2,0κ1,17+28κ2,02κ1,16+56κ2,03κ1,15+70κ2,04κ1,14+56κ2,05κ1,13+28κ2,06κ1,12+8κ2,07κ1,1+κ2,08 
κ1,18+28κ2,02κ1,16+70κ2,04κ1,14+28κ2,06κ1,12+κ2,08 κ1,19+9κ2,0κ1,18+36κ2,02κ1,17+84κ2,03κ1,16+126κ2,04κ1,15+126κ2,05κ1,14+84κ2,06κ1,13+36κ2,07κ1,12+9κ2,08κ1,1+κ2,09 

In the high-temperature approximation for the solvent fluctuations and for a simple Debye model of solvent relaxation, the correlation function of the energy gap becomes a single exponential.26 In this case, the double time integral with the energy gap correlation function in Eq. (19) can be approximated by a Gaussian,

0tdτ10τ1dτ2CE(τ1τ2)expλst2βħ2,
(22)

where λs is the solvent reorganization energy. Note that the simple Debye model of solvent relaxation breaks down when the short time scale solvent dynamics is important, in which case the energy gap correlation function may include Gaussian components to describe the short time scale relaxation processes.27,28 Using the readily available analytical expressions for the variance 〈δR2〉 and correlation function δR(0)δR(τ) of the quantum harmonic oscillator with mass M and frequency Ω,29 we obtain the following general expression for the reaction flux J(t):

J(t)=V02exp2λαħΩζexpn=11nλγζħΩn1nK̃n(1)(t)+λαζħΩK̃n+1(2)(t)×exp12χΩ2t2+p(cosΩt1)+i(qsinΩt+θΩt).
(23)

In this expression, λα=ħ2α22M and λγ=ħ2γ2M are the coupling reorganization energies corresponding to the linear and quadratic terms, respectively, in the logarithmic expansion of the coupling, ζ=cothβħΩ/2, and K̃n(i)=Kn(i)/δR2n for i = 1, 2. The remaining parameters are defined exactly as in our previous work:8 

(24)
p=ζλR+λαħΩ2λRλαħΩ,
q=λR+λαħΩ2ζλRλαħΩ,
χ=2λsβħ2Ω2,θ=ΔG+λsħΩ,
where λR=12MΩ2ΔR2 is the proton donor-acceptor mode reorganization energy, which is zero for systems with a symmetric proton transfer interface.

The additional exponential factor arising from the quadratic term in the coupling expansion is an infinite sign-alternating series that in certain regimes may cause convergence difficulties for the calculation of the reaction flux given in Eq. (23). The key temperature-dependent convergence parameter is ηγ = λγζ/ħΩ, which depends on γ as well as on the mass M and the frequency Ω of the proton donor-acceptor mode. In the low-frequency limit for the proton donor-acceptor mode, βħΩ ≪ 1, the convergence parameter becomes ηγ=γ/βMΩ2. Thus, at high temperatures and for a very soft proton donor-acceptor mode with small force constant MΩ2, the convergence may become slow. The convergence should improve significantly in the high-frequency limit, βħΩ ≫ 1, where ηγ=λγ/ħΩ1.

For the case in which only the linear term in the logarithmic coupling expansion is retained (i.e., γ = 0), closed analytical forms of the rate constant have been derived in the low- and high-frequency regimes for the proton donor-acceptor mode.5,8 In the case when the quadratic term is included (i.e., γ≠0), a closed analytical expression for the rate constant can be obtained if we neglect the equilibrium dynamics of the proton donor-acceptor coordinate R by replacing the time correlation function CR(t)=δR(0)δR(t) by its value at t = 0, CR(0)=δR2. This approximation for the proton donor-acceptor mode is valid when the characteristic time scale of the correlation function CR(t) is much longer than the time scale associated with the solvent damping term given in Eq. (22). In previous work,7 we analyzed the time dependence of the correlation function CR(t) calculated from molecular dynamics simulations of model PCET systems and its effects on the rate constant for the linear logarithmic coupling expansion. For typical PT interfaces, we found that the approximation CR(t) ≈ CR(0) is valid at ambient temperatures for solvents with large reorganization energy. In this case, the infinite series in Eq. (23) reduces to a time-independent prefactor, and the time integral of the reaction flux correlation function can be readily evaluated. For PT interfaces with Λ̃R=0, the rate constant acquires the following analytical form:

k=V02ħexp2λαζħΩ+2λγζ1+2λγζħΩ12πβλsexpβΔG+λs24λs.
(25)

The low-frequency (βħ Ω ≪ 1) and high-frequency (βħΩ ≫ 1) limits for the proton donor-acceptor mode can be obtained from the above expression by substituting the corresponding limits for ζ=cothβħΩ/2.

Unfortunately, in the general case, when the time dependence of CR(t) is retained and the interference effects between the short time dynamics of the proton donor-acceptor mode and the fluctuations of the energy gap are included, the inclusion of the quadratic term in the coupling expansion leads to expressions that are not easily simplified into this type of analytical form. In this case, the effects of the quadratic term on the rate constant cannot be represented as a simple prefactor in the rate constant expression. In our calculations of the nonadiabatic rate constants for model systems discussed below, we use Eq. (5) with the full reaction flux given in Eq. (23) and perform the integration numerically.

The alternative thermal averaging approach6 assumes that the proton donor-acceptor mode remains in thermal equilibrium during the reaction and that the nonadiabatic transitions or proton tunneling events occur independently at all values of R. This assumption is valid when the relaxation along R is fast compared to the rate of nonadiabatic transitions and is similar to the slow reaction limit in the Sumi-Marcus model for electron transfer.30 In this case, ignoring the weak dependence of the reorganization energy and the reaction free energy on R, the nonadiabatic rate constant is given by the following expres sion:

k=1ħV(R)2πβλsexpβΔG+λs24λs,
(26)

where the first prefactor is obtained by simple thermal averaging of the squared vibronic coupling over an equilibrium distribution of the proton donor-acceptor distances. For the exponential form of the coupling given in Eq. (4), the average is expressed as

V(R)2=P(R)V(R)2dR=V02P(R)exp2α(RR0)γ(RR0)2dR,
(27)

where P(R) is the equilibrium distribution function for proton donor-acceptor distances.

In the low-frequency (βħΩ ≪ 1) and high-frequency (βħΩ ≫ 1) limits for a harmonic oscillator representing the proton donor-acceptor vibrational mode, the integral in Eq. (27) can be evaluated analytically, resulting in the following analytical expressions:

V(R)2low-freq=V02exp2α2βMΩ2exp4α2γβMΩ2βMΩ2+2γ1+2γβMΩ21/2,
(28)
V(R)2high-freq=V02expħα2MΩexpħ2α2γMΩMΩ+ħγ1+ħγMΩ1/2.
(29)

As pointed out previously,31 in the absence of the quadratic term in the coupling expansion, the averaging of the squared coupling with the classical harmonic distribution function leads to a prefactor that coincides with the prefactor in the rate constant expression derived with the cumulant expansion approach in the low-frequency limit for the proton donor-acceptor mode. When the quadratic term is included in the coupling expansion, the expressions given in Eqs. (28) and (29) coincide with the low-frequency (βħΩ ≪ 1) and high-frequency (βħΩ ≫ 1) limits of the prefactor in Eq. (25). Thus, the thermal averaging of the rate constant with the equilibrium distribution function for the proton donor-acceptor mode, resulting in the rate constant expression in Eq. (26), corresponds to an approximation that completely neglects the effects from the dynamical interference between the proton donor-acceptor motion and the fluctuations of the energy gap. This relation between the cumulant expansion and thermally averaged rate constant expressions can be proven for any form of the vibronic coupling. From a theoretical perspective, the cumulant expansion approach is more general and consistent because the equilibrium dynamics along the R coordinate is treated explicitly on the same level as the dynamics of the energy gap.

The expression in Eq. (28) is similar to the result obtained by Trakhtenberg et al.1 for hydrogen tunneling reactions with a single promoting mode (cf. Eq. (8) in Ref. 1). Analyzing the tunneling splittings in proton potentials derived from Morse functions for chemical bonds, the authors concluded that the additional prefactors arising from the quadratic term in the logarithmic coupling expansion can be safely neglected. However, other studies3,12,32 have found that at high temperatures and for soft PT interfaces with very small force constants MΩ2, the second exponential factor in Eq. (28) can become small and thus can significantly reduce the rate constant. When only the linear term in the coupling expansion is included, the coupling is overestimated at short proton donor-acceptor distances, and for soft PT interfaces at high temperatures, the distribution function P(R) is broad enough to allow sampling of these short distances. On the other hand, the use of more realistic potentials along the coordinate R for hydrogen-bonded PT interfaces, such as Morse-like potentials with repulsive walls preventing the proton donor and acceptor from unphysical close contact, reduces the effects of the quadratic term due to decreased sampling probability at unphysical short donor-acceptor distances.

In an alternative, explicit thermal averaging approach, the vibronic coupling is calculated analytically or numerically at each proton donor-acceptor distance along a grid, followed by numerical evaluation of the first integral in Eq. (27).25,33–36 Thus, the vibronic couplings are calculated explicitly for the entire range of proton donor-acceptor distances sampled. This approach does not rely on the expansion of the vibronic coupling near the equilibrium proton donor-acceptor distance and therefore avoids any assumption of linear or quadratic behavior.

Another approach, which should be valid for a wide range of proton donor-acceptor mode effective masses M and frequencies Ω, employs a 2D quantum treatment of the PT interface. In this case, the R coordinate is treated on the same level as the electron and proton coordinates, and the reactant and product vibronic states are defined as solutions of the vibrational Schrödinger equation for the proton and its donor and acceptor moving on 2D diabatic electronic potentials UJ(q, R) corresponding to the reactant and product diabatic electronic states, J = I,  II. The shape of these diabatic potentials is assumed to be independent of the environmental coordinates ξ according to the definition of diabatic electronic states. In the case of 2D harmonic potentials UJ(q, R), the solutions of the vibrational Schrödinger equation can be obtained analytically,12 and the vibronic couplings can be evaluated analytically in the electronically nonadiabatic regime when they are proportional to the overlap integrals between the 2D vibrational wavefunctions, ψμJ(q,R), of the reactant and product states. A detailed analysis of such overlap integrals has been performed recently.12 For more realistic anharmonic 2D potentials, the 2D vibrational wavefunctions ψμJ(q,R) can be calculated numerically using, for example, Fourier Grid Hamiltonian (FGH) methods.37,38

Alternatively, the adiabatic separation between the proton coordinate q and the donor-acceptor coordinate R can be invoked by treating q as a fast coordinate and R as a slow coordinate. In this case, the reactant and product vibronic states are products of the adiabatic proton vibrational wavefunctions χμJ(q|R), parametrically dependent on the R coordinate, and vibrational wavefunctions for the R coordinate ϕμmJ(R). As mentioned in Section II A, such an approximation is justified for hydrogen-bonded PT interfaces due to the specific form of the coupling between the coordinates q and R and the large difference in associated masses.15 For adiabatically separated q and R in the reactant and product potentials, the nonadiabatic rate constant for a Debye solvent is given by

k=1ħμm{I}PμmIνn{II}ϕμmIVμν(R)ϕνnII2πβλsexpβΔGμm,νn+λs24λs,
(30)

where PμmI is the Boltzmann probability for the reactant state χμIϕμmI and ΔGμm,νn is the reaction free energy for the pair of vibronic states μm and νn.

In previous work,8 we compared the full 2D quantum treatment of coordinates q and R with other approaches for a series of models with anharmonic diabatic electronic potential energy surfaces. Numerical calculations of the nonadiabatic rate constants in the electronically nonadiabatic regime for various models illustrated that at ambient temperatures of ∼300 K and typical R-mode frequencies of 100–500 cm−1 and masses of 10–100 amu, the full 2D quantum and explicit thermal averaging approaches provide very similar results. Note that the general cumulant expansion approach described in Section II B2,5,8,10 also treats both the q and R coordinates quantum mechanically when the quantum mechanical correlation function for the harmonic R-mode is utilized. Moreover, similar to the thermal averaging approach, the 2D quantum approach neglects the interference effects arising from the dynamics of the proton donor-acceptor motion that are included in the cumulant expansion approach.

In this section, we analyze the effects of the quadratic term in the logarithmic expansion of the vibronic coupling on the rate constants for model systems spanning various regimes.

For the analysis in the electronically nonadiabatic regime, we employed two simple analytical models with harmonic and Morse proton potentials, denoted 1D-HARM and 1D-MORSE, respectively. These simple symmetric models assume that the shape of the diabatic proton potentials does not depend on the proton donor-acceptor coordinate R, which affects only the distance between the minima on the reactant and product potentials. The motion along R is assumed to be harmonic with effective mass M and frequency Ω. The harmonic proton potentials in the 1D-HARM model have an equilibrium distance of 1.0 Å, proton mass of 1.007 amu, and frequency of 2900 cm−1. The parameters of the Morse proton potentials in the 1D-MORSE model are similar to the parameters used to describe the nonadiabatic PCET reaction in soybean lypoxygenase:39 the bond dissociation energy D = 80 kcal/mol, the Morse parameter β = 2.1189 Å−1, and the bond equilibrium distance r0 = 1.0 Å. The proton vibrational wavefunctions and overlap integrals for the harmonic and Morse proton potentials were calculated analytically,40,41 and the coupling attenuation parameters α and γ were calculated by evaluating the derivatives in Eqs. (2) and (3) numerically.

For the analysis in the intermediate and electronically adiabatic regimes, we used 2D hybrid non-harmonic proton potentials that combine a London-Eyring-Polanyi-Sato (LEPS) potential ULEPS(q, R) for triatomic collinear systems and a harmonic potential Uharm(R) along the proton donor-acceptor coordinate R. The electronically diabatic potentials are represented by the diagonal elements of the 2 × 2 matrix of the electronic Hamiltonian with matrix elements:

HI,I(q,R)=UILEPS(q,R)+Uharm(R),HII,II(q,R)=UIILEPS(q,R)+Uharm(R),HI,II(q,R)=Vel,
(31)

where Vel is a constant electronic coupling parameter chosen to be 10 kcal/mol. The corresponding adiabatic potential U0(ad)(q,R) for the ground electronic state is obtained by diagonalization of this [2 × 2] matrix.

This potential was used in one of our earlier model studies,8 and the reader is referred to Ref. 8 for technical details and the values of the parameters for the LEPS potentials. The two symmetric models are denoted as LEPS-M10-Ω200 and LEPS-M50-Ω200 according to the values of the effective mass M and frequency Ω of the donor-acceptor mode. The contour plots for the diabatic 2D reactant and product potentials are shown in Figure 2 for one of the models. We point out that these 2D model potentials include the effects of Duschinski rotations that have been considered to be important.12 We used the FGH method for the numerical calculation of the tunneling splittings in the electronically adiabatic proton potential for one-dimensional slices along the proton coordinate q. From these tunneling splittings, we obtained the vibronic coupling in the electronically adiabatic limit and calculated the semiclassical vibronic couplings in the intermediate regime using Eq. (10). The corresponding coupling attenuation parameters α and γ were calculated numerically using Eqs. (2) and (3), respectively.

FIG. 2.

Contour plots of the reactant (left) and product (right) electronically diabatic potentials HI,I(q, R) and HII,II(q, R), respectively, for the LEPS-M50-Ω200 model. The effective mass and frequency of the proton donor-acceptor mode are M = 50 amu and Ω = 200 cm−1, respectively.

FIG. 2.

Contour plots of the reactant (left) and product (right) electronically diabatic potentials HI,I(q, R) and HII,II(q, R), respectively, for the LEPS-M50-Ω200 model. The effective mass and frequency of the proton donor-acceptor mode are M = 50 amu and Ω = 200 cm−1, respectively.

Close modal

In the electronically nonadiabatic regime, the vibronic coupling is of the form given in Eq. (15), which is proportional to the overlap integral between the reactant and product proton vibrational wavefunctions associated with the reactant and product diabatic electronic states, respectively. To estimate the significance of the quadratic term in the logarithmic coupling expansion, we calculated the linear and quadratic attenuation parameters, α and γ, respectively, for the proton vibrational overlap integrals for the model systems with harmonic and Morse proton potentials. The calculated values of these attenuation parameters for proton donor-acceptor distances ranging from 2.7 to 3.0 Å are provided in Table II. Because of the Gaussian nature of the harmonic overlap integrals, the γ attenuation parameters are larger for the harmonic potentials than for the Morse potentials. The asymmetry of the Morse potentials leads to slightly enhanced proton delocalization toward the acceptor site, leading to smaller values of the γ parameters.

TABLE II.

Attenuation parameters α (in Å−1) and γ (in Å−2) and coupling reorganization energies λα and λγ (in kcal/mol) for the overlap integrals between the ground state reactant and product vibrational wavefunctions for displaced harmonic and Morse proton potentials for hydrogen (H) and deuterium (D). The effective mass of the proton donor-acceptor mode used in the calculations of the coupling reorganization energies is M = 10 amu.

R0, Åα (H)γ (H)α (D)γ (D)λα (H)λγ (H)λα (D)λγ (D)
1D-HARM model 
2.7 30.32 43.32 42.87 61.25 4.432 0.209 8.860 0.295 
2.8 34.67 43.32 48.99 61.25 5.789 0.209 11.57 0.295 
2.9 38.99 43.32 55.12 61.25 7.327 0.209 14.65 0.295 
3.0 43.32 43.32 61.24 61.25 9.045 0.209 18.08 0.295 
1D-MORSE model 
2.7 19.83 20.64 28.69 29.18 0.379 0.0199 0.79 0.0281 
2.8 21.79 18.57 31.46 26.25 0.458 0.0179 0.95 0.0253 
2.9 23.55 16.70 33.95 23.61 0.535 0.0161 1.11 0.0228 
3.0 25.14 15.03 36.19 21.24 0.609 0.0145 1.26 0.0205 
R0, Åα (H)γ (H)α (D)γ (D)λα (H)λγ (H)λα (D)λγ (D)
1D-HARM model 
2.7 30.32 43.32 42.87 61.25 4.432 0.209 8.860 0.295 
2.8 34.67 43.32 48.99 61.25 5.789 0.209 11.57 0.295 
2.9 38.99 43.32 55.12 61.25 7.327 0.209 14.65 0.295 
3.0 43.32 43.32 61.24 61.25 9.045 0.209 18.08 0.295 
1D-MORSE model 
2.7 19.83 20.64 28.69 29.18 0.379 0.0199 0.79 0.0281 
2.8 21.79 18.57 31.46 26.25 0.458 0.0179 0.95 0.0253 
2.9 23.55 16.70 33.95 23.61 0.535 0.0161 1.11 0.0228 
3.0 25.14 15.03 36.19 21.24 0.609 0.0145 1.26 0.0205 

The important quantities affecting the nonadiabatic rate constant obtained by the cumulant expansion approach, based on the reaction flux given in Eq. (23), are the coupling reorganization energies λα and λγ, which are provided in Table II. The values of λγ are an order of magnitude smaller than the already very small values of λα for M = 10 amu, and they become even smaller as the effective mass of the proton donor-acceptor mode increases. The small magnitude of λγ results in small values of the convergence parameter ηγ that ensures convergence of the series in Eq. (23) at ambient temperatures. Figure 3(a) depicts the convergence parameter ηγ as a function of the effective mass M for R0 = 2.8 Å and Ω = 200 cm−1 at T = 300 K. The convergence of the series in the factor of the reaction flux that originates from the quadratic term in the coupling expansion requires only two or three terms for systems with moderate to large effective force constants associated with the proton donor-acceptor mode. For systems with small effective force constants (i.e., soft hydrogen-bonded PT interfaces), full convergence may require the inclusion of more terms. For example, to achieve convergence for the 1D-HARM model with M = 10 amu and Ω = 200 cm−1 at T = 300 K, we had to include terms up to tenth order in γ. For the more realistic Morse proton potentials, the convergence parameter is sufficiently small to ensure efficient convergence even for small effective force constants.

FIG. 3.

Key parameters in the rate constant expressions originating from the quadratic term in the vibronic coupling expansion, as functions of the effective mass M of the proton donor-acceptor vibrational mode. Solid and dashed lines are calculated using the overlap integrals of the ground state proton vibrational wavefunctions for Morse and harmonic reactant and product proton potentials, respectively. (a) The convergence parameter ηγ=λγcothβħΩ/2/ħΩ for the rate constant expression derived using the cumulant expansion approach. The frequency of the proton donor-acceptor mode used in these calculations is 200 cm−1. (b) The prefactors exp4α2γβMΩ2βMΩ2+2γ1+2γβMΩ21/2 and expħ2α2γMΩMΩ+ħγ1+ħγMΩ1/2 in the low-frequency (red lines) and high-frequency (blue lines) limits, respectively, for the rate constant obtained from the thermal averaging approach. The frequencies of the proton donor-acceptor mode used in these calculations are 200 cm−1 and 600 cm−1 for the low-frequency and high-frequency limits, respectively. The equilibrium proton donor-acceptor distance is R0 = 2.8 Å for all curves in both panels.

FIG. 3.

Key parameters in the rate constant expressions originating from the quadratic term in the vibronic coupling expansion, as functions of the effective mass M of the proton donor-acceptor vibrational mode. Solid and dashed lines are calculated using the overlap integrals of the ground state proton vibrational wavefunctions for Morse and harmonic reactant and product proton potentials, respectively. (a) The convergence parameter ηγ=λγcothβħΩ/2/ħΩ for the rate constant expression derived using the cumulant expansion approach. The frequency of the proton donor-acceptor mode used in these calculations is 200 cm−1. (b) The prefactors exp4α2γβMΩ2βMΩ2+2γ1+2γβMΩ21/2 and expħ2α2γMΩMΩ+ħγ1+ħγMΩ1/2 in the low-frequency (red lines) and high-frequency (blue lines) limits, respectively, for the rate constant obtained from the thermal averaging approach. The frequencies of the proton donor-acceptor mode used in these calculations are 200 cm−1 and 600 cm−1 for the low-frequency and high-frequency limits, respectively. The equilibrium proton donor-acceptor distance is R0 = 2.8 Å for all curves in both panels.

Close modal

For the thermally averaged rate constant, when an expansion of the coupling is used, the effects of the quadratic term on the rate constant can be significant for soft hydrogen-bonded PT interfaces with a small effective force constant. In this limit, the quadratic term is necessary to accurately describe the coupling at shorter proton donor-acceptor distances, which are sampled when the distribution function P(R) is broad (i.e., at high temperatures or for a very low-frequency proton donor-acceptor mode). As illustrated in Figure 3(b), however, the effects of the quadratic term become insignificant when the force constant of the proton donor-acceptor mode increases and the distribution function becomes more localized around the equilibrium value of the proton donor-acceptor distance in the reactant state. Note that the assumption of a harmonic potential along the proton donor-acceptor coordinate R is expected to become invalid for soft hydrogen-bonded interfaces. A more realistic anharmonic potential along the R coordinate would favor longer proton donor-acceptor distances, thereby avoiding the small distances at which the quadratic term in the coupling expansion is important.

For further analysis, we compared the rate constants obtained from the 2D quantum treatment of the proton donor-acceptor and proton coordinates, the thermally averaged rate constants with explicitly calculated couplings for all proton donor-acceptor distances, and the rate constants obtained from the cumulant expansion with and without the quadratic term. The results for M = 50 amu and M = 10 amu with Ω = 200 cm−1 are presented in Figure 4. For both models, the 2D quantum rate constants, thermally averaged rate constants with explicitly calculated couplings, and rate constants from the cumulant expansion with the quadratic term are in good agreement (dashed red, blue, and black lines, respectively). For the model with effective mass M = 50 amu, depicted in Figure 4(a), all of these rate constants are virtually indistinguishable for a wide temperature range. The rate constants calculated without the quadratic term in the coupling expansion are in good agreement with the others except for a slight deviation at higher temperatures (solid black line in Figure 4(a)), indicating that the effects of the quadratic term are insignificant for this model and therefore can be neglected.

FIG. 4.

Temperature dependence of the nonadiabatic rate constants for the proton donor-acceptor mode without the quadratic term (black line) and with the quadratic term (black dashed, blue dashed, and red dashed lines) for the 1D-MORSE model with R0 = 2.8 Å. The solid black line indicates the rate constant obtained from the cumulant expansion without the quadratic term (i.e., γ = 0). The dashed black lines indicate the rate constant obtained from the cumulant expansion including terms up to tenth order in γ. The red dashed lines indicate the rate constant obtained from the 2D quantum treatment of the proton donor-acceptor and proton coordinates. The blue dashed lines indicate the thermally averaged rate constant with the coupling calculated explicitly for the entire range of proton donor-acceptor distances. The dashed blue and red lines are virtually indistinguishable in both panels. The results are provided for models with (a) effective mass and frequency of the proton donor-acceptor mode M = 50 amu and Ω = 200 cm−1, respectively, and (b) effective mass and frequency of the proton donor-acceptor mode M = 10 amu and Ω = 200 cm−1, respectively.

FIG. 4.

Temperature dependence of the nonadiabatic rate constants for the proton donor-acceptor mode without the quadratic term (black line) and with the quadratic term (black dashed, blue dashed, and red dashed lines) for the 1D-MORSE model with R0 = 2.8 Å. The solid black line indicates the rate constant obtained from the cumulant expansion without the quadratic term (i.e., γ = 0). The dashed black lines indicate the rate constant obtained from the cumulant expansion including terms up to tenth order in γ. The red dashed lines indicate the rate constant obtained from the 2D quantum treatment of the proton donor-acceptor and proton coordinates. The blue dashed lines indicate the thermally averaged rate constant with the coupling calculated explicitly for the entire range of proton donor-acceptor distances. The dashed blue and red lines are virtually indistinguishable in both panels. The results are provided for models with (a) effective mass and frequency of the proton donor-acceptor mode M = 50 amu and Ω = 200 cm−1, respectively, and (b) effective mass and frequency of the proton donor-acceptor mode M = 10 amu and Ω = 200 cm−1, respectively.

Close modal

The model with a lower effective mass M = 10 amu, corresponding to a soft hydrogen-bonded PT interface, exhibits greater differences among the various rate constant expressions, as depicted in Figure 4(b). The deviation of the rate constant calculated from the cumulant expansion with the quadratic term (dashed black line in Figure 4(b)) from the 2D quantum and thermally averaged rate constants at higher temperatures is due to the slow convergence of the cumulant expansion. These results were generated by including terms up to tenth order in γ in the cumulant expansion rate constant, and the inclusion of higher-order terms is expected to improve the agreement at higher temperatures. On the other hand, the slight persistent deviation of the rate constant calculated from the cumulant expansion with the quadratic term (dashed black line in Figure 4(b)) from the 2D quantum and thermally averaged rate constants at ambient and lower temperatures is due to the neglect of the dynamical interference effects in the 2D quantum and thermally averaged rate constants. In this regime, the cumulant expansion rate constant is expected to be more accurate because it includes the dynamical interference between the proton donor-acceptor mode and the energy gap.

More significant differences are observed between the rate constants calculated without the quadratic term in the coupling expansion (solid black line in Figure 4(b)) and the rate constants that include the quadratic term, implying that the quadratic term is important when shorter proton donor-acceptor distances are more accessible. Note that in realistic chemical or biological systems, the shorter proton donor-acceptor distances will not be as readily accessible even for weakly hydrogen-bonded PT interfaces due to strong repulsive nonbonding interactions between the donor and acceptor groups. Thus, the models based on a harmonic description of a low-frequency proton donor-acceptor mode appear to overemphasize the region of short distances that are not relevant for physically realistic systems.

In the electronically adiabatic regime, typical for PT and HAT reactions, the vibronic coupling V(R) at each proton donor-acceptor distance R is related to the tunneling splitting for a symmetric double well proton potential on the electronically adiabatic ground state surface U0(ad)(q,R). Because the models studied herein are symmetric, the proton potentials are symmetric for a fixed proton donor-acceptor distance R. Calculation of the semiclassical couplings and associated electron-proton adiabaticity parameters for a range of proton donor-acceptor distances for the two models LEPS-M10-Ω200 and LEPS-M50-Ω200 with different effective masses M = 10 amu and M = 50 amu, respectively, revealed that both models are in the intermediate regime. This property is illustrated in Figure 5 for the LEPS-M50-Ω200 model.

FIG. 5.

Proton donor-acceptor distance dependence of the vibronic coupling for the LEPS-M50-Ω200 model calculated using the electronically nonadiabatic limit expression (blue dashed line), the electronically adiabatic limit expression (red dashed line), and the semiclassical expression bridging both limits (black line).

FIG. 5.

Proton donor-acceptor distance dependence of the vibronic coupling for the LEPS-M50-Ω200 model calculated using the electronically nonadiabatic limit expression (blue dashed line), the electronically adiabatic limit expression (red dashed line), and the semiclassical expression bridging both limits (black line).

Close modal

The equilibrium proton donor-acceptor distances R0 and attenuation parameters α and γ for the overlaps, tunneling splittings, and semiclassical vibronic couplings are given in Table III. The values of the attenuation parameters in these models are similar to those calculated for the overlaps in the 1D-MORSE model (Table II). Thus, the effects of the quadratic term in the coupling expansion on the rate constants in the various limits are also expected to be similar.

TABLE III.

Attenuation parameters α (in Å−1) and γ (in Å−2) for the overlap integrals, tunneling splittings, and semiclassical vibronic couplings for the two LEPS models.

OverlapTunneling splittingSemiclassical couplinga
R0, Åαγαγαγ
LEPS-M10-Ω200 model 
2.72 17.81 25.43 16.77 27.67 16.87 27.23 
LEPS-M50-Ω200 model 
2.78 19.75 25.84 19.33 26.45 19.39 26.08 
OverlapTunneling splittingSemiclassical couplinga
R0, Åαγαγαγ
LEPS-M10-Ω200 model 
2.72 17.81 25.43 16.77 27.67 16.87 27.23 
LEPS-M50-Ω200 model 
2.78 19.75 25.84 19.33 26.45 19.39 26.08 
a

The effective proton tunneling time is τp = 0.4479 (0.4267) fs, the effective electronic transition time is τe = 1.518 (1.518) fs, and the electron-proton adiabaticity parameter is p = 0.2951 (0.2811) for the LEPS-M10-Ω200 (LEPS-M50-Ω200) models.

To illustrate the effects on the thermally averaged rate constant, we examined the distance dependence of the semiclassical vibronic coupling and the approximations to this coupling based on the logarithmic expansion, as given in Eq. (4), in relation to the width of the classical distribution function for the proton donor-acceptor distances. The classical distribution function P(R)expε0I(R)/kBT was calculated for the ground state reactant vibronic potentials ε0I(R) obtained by quantization of the proton moving on the electronically diabatic reactant electronic surface HI,I(q, R). The results are depicted in Figure 6 for two models with soft and rigid hydrogen-bonded PT interfaces. In both cases, the second-order logarithmic expansion of the coupling is sufficient to reproduce the semiclassical vibronic coupling nearly quantitatively. For the soft R-mode (Figure 6(b)), significant deviations of the first-order approximation to the coupling from the semiclassical vibronic coupling are observed within the sampling region covered by the distribution function P(R). In contrast, for the more rigid R-mode (Figure 6(a)), the first-order approximation to the vibronic coupling is satisfactory within the sampling region covered by the distribution function. Thus, the quadratic term will impact the rate constant only for the case of PT interfaces with weak hydrogen bonds and soft proton donor-acceptor vibrational modes.

FIG. 6.

Proton donor-acceptor distance dependences of the semiclassical vibronic coupling (black line) and the approximations to this coupling based on the logarithmic expansion to first order (blue dashed line) and second order (red dashed line), as given in Eq. (4). The results were obtained for (a) the LEPS-M50-Ω200 model and (b) the LEPS-M10-Ω200 model. The gray shaded Gaussian is the classical distribution function (in arbitrary units) for the proton donor-acceptor distances associated with the ground reactant vibronic state for each of the models.

FIG. 6.

Proton donor-acceptor distance dependences of the semiclassical vibronic coupling (black line) and the approximations to this coupling based on the logarithmic expansion to first order (blue dashed line) and second order (red dashed line), as given in Eq. (4). The results were obtained for (a) the LEPS-M50-Ω200 model and (b) the LEPS-M10-Ω200 model. The gray shaded Gaussian is the classical distribution function (in arbitrary units) for the proton donor-acceptor distances associated with the ground reactant vibronic state for each of the models.

Close modal

As discussed above, more realistic potentials that include repulsion between the proton donor and acceptor would prevent the sampling of the unphysical short distances. In this case, the short distances sampled extensively for the soft R-mode system described by Figure 6(b) would be energetically inaccessible, and the probability distribution P(R) would not be a Gaussian function but rather would be steeper at short distances. As in the nonadiabatic case, the thermal averaging procedure can also be performed with the explicit calculation of the vibronic coupling for each proton donor-acceptor distance followed by numerical integration, thereby avoiding the expansion of the coupling altogether and allowing a non-Gaussian probability distribution function.

We also emphasize that the Golden Rule formalism, which is based on second-order perturbation theory, is invalid for short proton donor-acceptor distances associated with large vibronic coupling. The Golden Rule formalism is the basis for all of the rate constant expressions discussed in this paper, including those in the electronically adiabatic and nonadiabatic limits, as well as the intermediate regime. Thus, all of these nonadiabatic rate constants are expected to be invalid when the thermal averaging procedure includes significant contributions from such short proton donor-acceptor distances. Models that allow the sampling of unphysical short distances should be treated with great care within these frameworks.

In this paper, we presented rate constant expressions for vibronically nonadiabatic PT and PCET reactions using a cumulant expansion approach based on the Golden Rule formalism. The resulting rate constant expressions are valid in both the electronically adiabatic and electronically nonadiabatic regimes, as well as the intermediate regime, as long as the appropriate form of the vibronic coupling is utilized. The general semiclassical expression for the vibronic coupling and its form in the electronically adiabatic and nonadiabatic limits are available from the literature. These rate constants differ from the previous rate constants derived with the cumulant expansion approach in that the logarithm of the vibronic coupling was expanded to second order. Thus, the logarithmic expansion of the vibronic coupling includes a quadratic as well as a linear term. The significance of this quadratic term was explored for a series of model systems.

Our analysis illustrates that inclusion of the quadratic term in the logarithmic expansion of the nonadiabatic coupling may lead to significant effects on the PT or PCET rate constant at high temperatures for PT interfaces with weak hydrogen bonds and soft proton donor-acceptor vibrational modes. The inclusion of the quadratic term in the framework of the cumulant expansion form of the nonadiabatic rate constant leads to a sign-alternating series in the expression for the time-dependent reaction flux correlation function. The convergence of this series may become slow at high temperatures for systems with a small effective force constant associated with the proton donor-acceptor mode. The effects of the quadratic term may also become significant when using the thermal averaging procedure for calculating the nonadiabatic rate constants in these regimes.

In both the cumulant expansion and the thermally averaged formulations, the effects of the quadratic term for weak hydrogen-bonding systems are less significant for more physically realistic models that prevent the sampling of unphysical short proton donor-acceptor distances. In particular, the assumption of a harmonic potential along the proton donor-acceptor coordinate is invalid for soft hydrogen-bonded interfaces, and more realistic anharmonic potentials along this coordinate inhibit sampling of the small proton donor-acceptor distances at which the quadratic term is important. Moreover, all of these rate constant expressions are based on the Golden Rule formalism, which is invalid at short proton donor-acceptor distances associated with large vibronic couplings. Thus, these rate constant expressions should be used only for systems that remain in the vibronically nonadiabatic limit (i.e., the Golden Rule regime) at all proton donor-acceptor distances sampled significantly. According to our analysis, the inclusion of the quadratic term in the coupling expansion is less likely to substantially impact the PT or PCET rate constant for physically realistic systems that remain in the Golden Rule regime.

In addition to the analysis of the quadratic term, the rigorous relation between the cumulant expansion and thermal averaging approaches was clarified. The cumulant expansion rate constant was shown to include effects arising from dynamical interference between the proton donor-acceptor and solvent motions. These dynamical effects are not included in the thermal averaging approach or in the two-dimensional full quantum approach. The general cumulant expansion and thermal averaging rate constants become mathematically equivalent for any form of the vibronic coupling when these dynamical effects are neglected. Understanding the regimes in which each rate constant expression is valid is important for future applications to PT and PCET in chemical and biological processes.

We are grateful for support from National Institutes of Health Grant No. GM056207 (applications to enzymes) and the Center for Molecular Electrocatalysis, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (applications to molecular electrocatalysts).

1.
L. I.
Trakhtenberg
,
V. L.
Klochikhin
, and
S. Y.
Pshezhetsky
,
Chem. Phys.
69
,
121
(
1982
).
2.
D.
Borgis
,
S.
Lee
, and
J. T.
Hynes
,
Chem. Phys. Lett.
162
,
19
(
1989
).
3.
Z. K.
Smedarchina
,
Chem. Phys.
150
,
47
(
1991
).
4.
A.
Suarez
and
R. J.
Silbey
,
J. Chem. Phys.
94
,
4809
(
1991
).
5.
D.
Borgis
and
J. T.
Hynes
,
Chem. Phys.
170
,
315
(
1993
).
6.
A. M.
Kuznetsov
and
J.
Ulstrup
,
Can. J. Chem.
77
,
1085
(
1999
).
7.
E.
Hatcher
,
A.
Soudackov
, and
S.
Hammes-Schiffer
,
Chem. Phys.
319
,
93
(
2005
).
8.
A.
Soudackov
,
E.
Hatcher
, and
S.
Hammes-Schiffer
,
J. Chem. Phys.
122
,
014505
(
2005
).
9.
M. V.
Basilevsky
and
V. A.
Tikhomirov
,
Mol. Phys.
106
,
2391
(
2008
).
10.
D.
Borgis
and
J. T.
Hynes
,
J. Phys. Chem.
100
,
1118
(
1996
).
11.
L.
Lavtchieva
and
Z.
Smedarchina
,
Chem. Phys. Lett.
184
,
545
(
1991
).
12.
A.
Benabbas
,
B.
Salna
,
J. T.
Sage
, and
P. M.
Champion
,
J. Chem. Phys.
142
,
114101
(
2015
).
13.
R. R.
Dogonadze
,
A. M.
Kuznetsov
, and
M. A.
Vorotyntsev
,
Phys. Status Solidi B
54
,
125
(
1972
).
14.
R. R.
Dogonadze
,
A. M.
Kuznetsov
, and
M. A.
Vorotyntsev
,
Phys. Status Solidi B
54
,
425
(
1972
).
15.
M. V.
Vener
and
N. D.
Sokolov
,
Chem. Phys. Lett.
264
,
429
(
1997
).
16.
Y.
Georgievskii
and
A. A.
Stuchebrukhov
,
J. Chem. Phys.
113
,
10438
(
2000
).
17.
M. D.
Newton
,
Chem. Rev.
91
,
767
(
1991
).
18.
M. S.
Child
,
Molecular Collision Theory
(
Academic Press
,
New York
,
1974
).
19.
B.
Fain
,
Theory of Rate Processes in Condensed Media
(
Springer
,
New York
,
1980
).
20.
S.
Hammes-Schiffer
,
Energy Environ. Sci.
5
,
7696
(
2012
).
21.
J. P.
Layfield
and
S.
Hammes-Schiffer
,
Chem. Rev.
114
,
3466
(
2014
).
22.
P. M.
Kiefer
and
J. T.
Hynes
,
J. Phys. Chem. A
107
,
9022
(
2003
).
23.
P. M.
Kiefer
and
J. T.
Hynes
,
J. Phys. Chem. A
108
,
11793
(
2004
).
24.
P. M.
Kiefer
and
J. T.
Hynes
,
J. Phys. Chem. A
108
,
11809
(
2004
).
25.
B.
Auer
,
L. E.
Fernandez
, and
S.
Hammes-Schiffer
,
J. Am. Chem. Soc.
133
,
8282
(
2011
).
26.
L. D.
Zusman
,
Chem. Phys.
49
,
295
(
1980
).
27.
J. T.
Hynes
,
J. Phys. Chem.
90
,
3701
(
1986
).
28.
J. T.
Hynes
and
M.
Bruehl
,
J. Phys. Chem.
96
,
4068
(
1992
).
29.
D. W.
Chandler
,
Introduction to Modern Statistical Mechanics
(
Oxford University Press
,
New York, Oxford
,
1987
), Vol.
1
.
30.
H.
Sumi
and
R. A.
Marcus
,
J. Chem. Phys.
84
,
4894
(
1986
).
31.
S.
Hammes-Schiffer
,
E.
Hatcher
,
H.
Ishikita
,
J. H.
Skone
, and
A. V.
Soudackov
,
Coord. Chem. Rev.
252
,
384
(
2008
).
32.
Z. K.
Smedarchina
and
W.
Siebrand
,
Chem. Phys.
170
,
347
(
1993
).
33.
M. P.
Meyer
and
J. P.
Klinman
,
Chem. Phys.
319
,
283
(
2005
).
34.
A.
Kohen et al.
, in
Hydrogen-Transfer Reactions
, edited by
J. T.
Hynes et al.,
(
WILEY-VCH Verlag GmbH & Co.
,
Weinheim
,
2007
), p.
1311
.
35.
S. J.
Edwards
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Phys. Chem. A
113
,
2117
(
2009
).
36.
D.
Roston
,
C. M.
Cheatum
, and
A.
Kohen
,
Biochemistry
51
,
6860
(
2012
).
37.
C. C.
Marston
and
G. G.
Balint-Kurti
,
J. Chem. Phys.
91
,
3571
(
1989
).
38.
G. G.
Balint-Kurti
and
C. L.
Ward
,
Comput. Phys. Commun.
67
,
285
(
1991
).
39.
E.
Hatcher
,
A. V.
Soudackov
, and
S.
Hammes-Schiffer
,
J. Am. Chem. Soc.
129
,
187
(
2007
).
40.
J. P.
Dahl
and
M.
Springborg
,
J. Chem. Phys.
88
,
4535
(
1988
).
41.
F.
Iachello
and
M.
Ibrahim
,
J. Phys. Chem. A
102
,
9427
(
1998
).