Nonadiabatic quantum transition-state theory in the golden-rule limit. II. Overcoming the pitfalls of the saddle-point and semiclassical approximations

We describe a path-integral molecular dynamics implementation of our recently developed golden-rule quantum transition-state theory (GR-QTST). The method is applied to compute the reaction rate in various models of electron transfer and benchmarked against exact results. We demonstrate that for systems exhibiting two or more transition states, rates computed using Wolynes theory [P. G. Wolynes, J.\ Chem.\ Phys.\ 87, 6559 (1987)] can be overestimated by orders of magnitude, whereas the GR-QTST predictions are numerically accurate. This is the case both at low temperature, where nuclear tunneling makes a considerable contribution, and also in the classical limit, where only GR-QTST rigorously tends to the correct result. Analysis shows that the saddle-point approximation employed by Wolynes theory is not valid in this case, which results in predictions of unphysical reaction pathways, whilst the energy constraint employed by GR-QTST resolves this problem. The GR-QTST method is also seen to give accurate results for a strongly anharmonic system by sampling configurations around the instanton pathway without making the semiclassical approximation. These promising results indicate that the GR-QTST method could be an efficient and accurate approach for simulating electron-transfer reactions in complex molecular systems.


I. INTRODUCTION
2][3][4] Such reactions can be described by a system of two weakly-coupled electronic states interacting with a large number of nuclear degrees of freedom. 5Electrontransfer reactions typically take place in the golden-rule limit and one cannot therefore employ the Born-Oppenheimer approximation.Instead, the rate is formally given by Fermi's golden-rule. 6In practice, however, further approximations are also required for a system in which the exact nuclear wave functions cannot be calculated.Marcus theory, 7,8 which effectively treats the free-energy curves harmonically and neglects nuclear quantum effects, has been widely applied to model these reactions.[30][31][32] We will concentrate on methods based on the path-integral description of quantum mechanics 33 which are able to describe nuclear tunneling effects.5][36] Both theories can be derived in a similar way using a number of steepest-descent approximations to the path-integral formulation of the flux correlation function to obtain a formula for the rate defined in terms of imaginary-time periodic classical trajectories and harmonic fluctuations around them. 18,37,38Each of these periodic trajectories is called an "instanton", and together they define a set of optimal tunneling pathways, each of which is identified as a mechanism of the reaction.Like the original version, [39][40][41] the golden-rule instanton theory can also be efficiently evaluated in the ring-polymer formulation, 19 and has been benchmarked on system-bath models and shown to give accurate predictions for electron-transfer rates. 20wever, although the theory is applicable to anharmonic systems, as the semiclassical limit only treats fluctuations to second order, it is equivalent to performing a local harmonic approximation, which cannot be used successfully in atomistic solvent environments.
Wolynes theory, 9 on the other hand, does not make a semiclassical approximation and instead evaluates the path integral numerically by sampling discretized paths, known as ring polymers. 42It can therefore be applied to atomistic environments using techniques developed for classical molecular dynamics. 43,44However, only imaginary-time (statistical) properties are rigorously available from this approach and so an expression for the rate is derived by integrating out the real-time behaviour using a saddle-point approximation.The method has been employed to compute electron-transfer rates in a number of applications. 28,45 Ref. 46 (hence referred to as Paper I), we proposed a new golden-rule quantum transition-state theory (GR-QTST) for predicting rates of electron-transfer reactions.This method was designed to go beyond instanton theory by sampling paths around the instanton in a similar way to ring-polymer transition-state theory 40,47,48 and the projected quantum instanton method. 49It also extends the ideas of classical golden-rule transition-state theory, 5,50 in which energy conservation from the reactant to product states is enforced by a constraint.
The method thus differs from Wolynes theory in that a constraint is applied to every path in the ensemble.In this respect, there is some similarity to the kinetically-constrained ring-polymer molecular dynamics (KC-RPMD) method, 16,17,51 although there are important differences in the definition and implementation.
In the previous work, 46 we tested the GR-QTST method on a spin-boson model and a one-dimensional anharmonic system and found that the results were consistently in good agreement with the exact rate and were also more accurate than published KC-RPMD results. 16,17In addition, we tested other approaches including instanton theory and Wolynes theory, which were also seen to perform well.In fact, Wolynes and instanton theories typically gave slightly more accurate results than GR-QTST in these cases (except for a particular discrepancy of Wolynes theory, which does not tend to the correct classical limit for high temperatures in an anharmonic system).However, we will show in Sec.IV that these two methods break down when applied to more complex systems than the simple models employed in Paper I, leaving GR-QTST as the only consistently reliable method.
In this work, we argue that the saddle-point approximation of Wolynes theory is not an asymptotic approximation in , and for this reason does not rigorously tend to the correct classical limit.We present a simple example where it completely breaks down by introducing a system with two transition states in Sec.II.Here, transition states (TSs) are defined loosely from a classical viewpoint as a point of minimum energy on the crossing seam between two diabatic surfaces. 50We note that Wolynes himself already hinted at the possibility of a break-down in such systems. 9By "break down" we mean that the method overestimates or underestimates the rate by at least an order of magnitude.A rate method liable to make this sort of error will also typically lose the ability to correctly identify the reaction mechanism, and so we will also evaluate the utility of the various methods with respect to this important criterion.
In order to perform the GR-QTST calculations presented in Sec.IV, we introduce a path-integral molecular dynamics (PIMD) implementation in Sec.III.For the strongly anharmonic systems employed in this work, this approach is more efficient than the importance sampling method employed in Paper I and will in principle allow the method to be applied to atomistic simulations.
It is clear that the semiclassical approximation of instanton theory breaks down when the fluctuations around the instanton pathway become strongly anharmonic, such as in a system with explicit solvents.It was exactly for this reason that we developed the GR-QTST method, and we demonstrate in Sec.IV D that the new approach remains accurate even when instanton theory fails.Finally, we conclude in Sec.V that GR-QTST successfully overcomes the pitfalls of both the saddle-point and semiclassical approximations.

II. THEORY AND ANALYSIS
In this section, we discuss semiclassical instanton theory and Wolynes theory, two methods derived in different ways from saddle-point approximations.The limitations of instanton theory due to its harmonic treatment of path fluctuations are well understood, but Wolynes theory has not been subjected to the same scrutiny.In particular, we will explain how the saddle-point approximation employed in the derivation of Wolynes theory can lead to errors greater than an order of magnitude in the rate predictions, even if the system is in the classical limit.We then show how the energy constraint introduced in the GR-QTST method corrects for the problems of both methods such that accurate results are obtained.
Consider an electron-transfer process from the reactant diabatic state |0 to the product diabatic state |1 . 5The Hamiltonian is where ∆ is the electronic coupling, and Ĥn = p2 /2m + V n (x) is a one-dimensional nuclear Hamiltonian for the electronic state |n with potential-energy surface V n (x).The nuclear coordinate is x which has associated mass m and conjugate momentum p.Because our methods are based on a path-integral formalism, they are easily extended to treat a system with D nuclear degrees of freedom with nuclear coordinates x = (x 1 , . . ., x D ) and conjugate momenta p = (p 1 , . . ., p D ), each with an associated mass m j .For these multidimensional systems, the nuclear Hamiltonian for each diabatic state is Ĥn = D j=1 p2 j /2m j +V n (x).Here we will assume that ∆ is a constant, which is known as the Condon approximation, but all the methods we present are easily generalized to allow for a position-dependent coupling.
3][54] We shall assume that the reaction proceeds in the nonadiabatic golden-rule limit, i.e. ∆ is small enough such that only the leading term need be considered for the rate. 55The correlation function then takes on a particularly simple form and leads to an expression for Fermi's golden-rule rate constant given by 57 where The symbol Tr n is a trace over nuclear degrees of freedom only.Typically this trace cannot be evaluated exactly.However, in this section, we will illustrate the theoretical concepts using a simple model system where the analytical result can be found (Appendix A).
It should be noted that this exact expression for the rate is independent of the choice of τ .This is because φ u (z) is an analytic function of the complex variable z, where τ = Re z and t = − Im z.Therefore the contour integral can be deformed to effectively change the value of τ without affecting the result of the integral.Finally, note that φ u (z) is a real number when t = 0, but is complex in general.
In general, the major difficulty in evaluating the integral in Eq. ( 2) numerically by quadrature is because the integrand itself (and its first few derivatives) can only be obtained at short times, as imaginary-time path-integral sampling techniques only allow an efficient numerical computation for t = 0. 42 We will discuss two methods in which this integral is carried out by the method of steepest descent, which is also known as a saddle-point approximation. 58e semiclassical instanton rate 18 is obtained by replacing the two propagators in Eq. ( 3) by their semiclassical limits (a sum over imaginary-time classical paths) 59 and evaluating the trace over nuclear degrees of freedom by steepest-descent integration.We call this the semiclassical approximation because it is exact in the limit → 0 (keeping the total imaginary time, β , fixed).Finally a further steepest-descent approximation is taken for the time integral.Alternatively, one can perform the steepest-descent integrals over path variables and time simultaneously, leading to the same result. 38The stationary points in this extended space correspond to instanton configurations, each with t = 0 but a different value of τ .Taking the high-temperature, large-mass limit gives a formula 18 which is recognized as a local harmonic approximation to classical golden-rule transition-state theory. 5,50lynes theory 9 does not employ the semiclassical approximation (steepest-descent integration over the path variables) to compute the effective action, but does employ a saddlepoint approximation for the integral over time in Eq. ( 2).This approximation is equivalent to expanding φ u (z) as a Taylor series to second order about its saddle point, thus effectively assuming that the correlation function has a Gaussian form centred on t = 0, whose integral is known.This is commonly called a second-order cumulant expansion.The resulting expression for the rate is formulated as where τ * is a real number and is defined as the stationary point which maximizes the unconstrained effective action, φ u (τ ). 60This approximation has been called a saddle-point approximation, 9 but note that is not an asymptotic approximation 58 in as φ u (τ − it) itself depends on .This is in contrast to the steepest-descent integration used to derive semiclassical instanton theory, which is a rigorous asymptotic approximation 37,38,61 and thus becomes exact in the limit → 0.
As pointed out in Ref. 50 and demonstrated numerically in Paper I, the Wolynes rate does not tend correctly to the classical limit at very high temperatures (unless the system is harmonic).Nonetheless, the errors in previous examples were rather small.In this paper, we demonstrate that orders-of-magnitude error can result from its saddle-point approximation, in particular for systems with more than one transition state.
We shall investigate a system which exhibits two diabatic crossings, denoted as A and B.
If we assume that the two crossing points are well separated in coordinate space, or couple different electronic states, the effective action is e −φu(τ +it)/ = e −φ (A) u (τ +it)/ + e −φ (B) u (τ +it)/ . ( This expression is analogous to forming a joint probability distribution by algebraically adding two independent probability distributions.The exact rate expression can therefore be broken into two parts, k = k (A) + k (B) , where k (s) is the rate through one of the crossings s ∈ {A, B}.Because the integration of a sum is simply the sum of integrals, it can be shown that the classical and semiclassical instanton rates are rigorously separable in this same way.
It is also clear that this same argument can be generalized to systems with more than two crossings.
In contrast, Wolynes theory is not separable and the rate obtained from one single calculation is different from the sum of two separate calculations.In fact, if we apply Wolynes theory directly to the total effective action, φ u (τ ), the saddle-point approximation may be of very poor quality.This is in spite of the fact that the sum of two separate Wolynes theory calculations on the individual crossings may give a much better approximation to the rate.
We call this modified method "separated Wolynes theory".curves are evaluated separately using Eq.(A2) for the transition states corresponding to A (blue) and B (green), whereas φ u (τ ) (dashed red) is the combined result, defined by Eq. ( 5).
We demonstrate these ideas by considering the specific example of a system with a harmonic reactant potential weakly coupled to two product states described by linear potentials which make crossings at the same height but with different slopes.The parameters are defined in Appendix A along with analytic expressions for the unconstrained effective action.u (τ ) which are compared to φ u (τ ), defined by Eq. ( 5) with t = 0.It is seen that each of these three action functions has a maximum at a different value of τ .
Therefore, following the original Wolynes theory, 9 one would obtain a single value of τ * at the maximum of φ u (τ ) and employ a saddle-point approximation around this value to yield the rate.However, this does not lead to a good approximation because, as is shown in Fig. 2(a), the correlation function for this system defined with τ = τ * can be strongly oscillatory.In principle, the exact rate can be obtained as the integral this correlation function, but the cumulant expansion employed by Wolynes theory attempts to approximate it by a Gaussian.The naive saddle-point approximation is therefore clearly not valid in the usual sense of a steepest-descent approximation, which is supposed to remove the imaginary component of the exponent, 58 but here does so only for a very small region around the origin.
Wolynes theory, which can only approximate the positive part of the first peak, can thus overestimate the integral by many orders of magnitude.
The correlation functions for the two transition states can also be investigated separately.
The exact total rate is given by the sum of the integrals over one blue and one green curve chosen from either Fig. 2(b) or (c).One can see that either of the e −φ (s) u (τ +it)/ correlation functions can be made approximately Gaussian, but that different choices of τ * s are required for this, which cannot be simultaneously satisfied.The choice τ = τ * made by Wolynes theory removes neither the oscillations of e −φ (A) u (τ +it) nor of e −φ (B) u (τ +it) .The semiclassical approximation of instanton theory follows a different approach and integrates over the blue part of (b) and the green of (c), both of which can be well approximated by Gaussians.In principle Wolynes theory could be systematically improved by going to higher-order cumulant expansions, but this would not be practical because a huge number of terms would be required to describe the oscillatory behaviour.It is therefore clear that in this case, the semiclassical approximation is superior to the cumulant expansion.However, note that this may not be the case for strongly anharmonic systems such as those we consider in Sec.IV.
A physical interpretation for the derivative of the effective action is a measure of the difference in energy between reactant and product states.Thus choosing τ = τ * s at the stationary point of φ (s) u (τ ) ensures that the energies match around this transition state, as is required by Fermi's golden rule. 50However, the value of τ = τ * which maximizes φ u (τ ) matches the average reactant energy to the average product energy globally, but they are not necessarily matched locally around each transition state.Wolynes theory will therefore sample unphysical paths that do not resemble instanton-like configurations, leading to a loss of mechanistic insight as well as orders-of-magnitude error in the predicted rate.We can also make some statements about the effect of this problem in more general systems.Firstly, assuming that both φ (s) u curves exhibit only one maximum each and that τ * A < τ * B , then τ * , which is the maximum of φ u (τ ), will always be in between the two maxima, i.e. τ * A < τ * < τ * B . 62Also, φ u (τ * ) will always be smaller than either of φ (s) u (τ * s ).Therefore, although Wolynes theory is not variational in general, the problem we identify here will typically cause the rate to be overpredicted.However, it will only have a significant impact in cases for which φ A and τ * B .Of course, for this simple model system, one can fix the problem by computing the rates using the separated Wolynes method, in which the saddle-point approximation is evaluated independently for the two transition states at τ * A and τ * B .The total rate is given by the sum of these two rates.Unfortunately it is not always so easy to separate the different contributions to Wolynes theory, in particular for complex molecular environments.It is for this reason we have proposed the GR-QTST method, 46 which we have developed to avoid this problem.
The results from various rate calculations are presented in Table I, where they can be compared with the benchmark exact result from Fermi's golden-rule.The difference between classical and exact rates for this system indicates that we are in the quantum regime where tunneling is responsible for enhancing the rates by over three orders of magnitude.It is also worth noting that due to its neglect of tunneling, classical mechanics incorrectly predicts a faster rate through transition state A. Semiclassical instanton rates, however, are very accurate, which is expected for this system which includes no anharmonicities such that the semiclassical approximation is exact for integrals over nuclear degrees of freedom.Therefore, there is only a very minor difference between semiclassical instanton theory and the separated Wolynes theory here.In both cases the only approximation used is to employ the saddlepoint method to the time integral over each individual Gaussian-like correlation function.
Instanton theory is calculated separately at τs , which is the maximum of the classical action S (s) (τ ) [Eq. (A3)], and separated Wolynes theory at τ * s , which is the maximum of φ (s) u (τ ) [Eq. (A2)].For this system, we find that τs and τ * s differ by less than 1% at either of the transition states, thus giving almost identical predictions from the two theories.
Despite the success of the separated version, standard Wolynes theory clearly overestimates the rates by over two orders of magnitude.In fact, this problem persists even in the classical limit.To demonstrate this, we calculated rates on the 1D system over a range of "quantumness", defined by parameter α, as shown in Fig. 3. Small α values correspond to the classical limit, while for large α values nuclear quantum effects are prominent.We find that the standard Wolynes rate, which is estimated by treating both TSs together, incorrectly predicts a huge speed-up in the rate relative to a classical calculation, even in the α → 0 limit.The error in log k is actually largest in the classical limit for this system because the two TSs in the model make almost equal contributions to the classical rate, GR-QTST [Eq.( 13)] and Wolynes calculations [Eq.( 4)] for the one-dimensional system defined by Eq. (A1).N = 200 ring-polymer beads were used to perform the GR-QTST calculations with 5000 Monte-Carlo steps.The statistical errors of one standard deviation are given.Instanton and Wolynes rate calculations have been performed analytically using φ Comparison between an instanton calculation with a finite number of beads and the analytic result shows that N = 200 is sufficient to describe the tunneling accurately.To make the results dimensionless, all the rates have been divided by k cl , which is the total classical rate for the reaction through both transition states.which causes the φ u (τ ) curve to deviate most strongly from either of the separated effective actions in this limit.Note however that in the quantum regime, even though the exact rate is completely dominated by transition state A, the action φ u (τ ) is still significantly perturbed by the presence of the other TS, leading to a large error in Wolynes theory.
GR-QTST is not based on a saddle-point approximation, but on its connection to instanton theory.The main difference to Wolynes theory is that the rate is obtained from a constrained effective action, φ c (τ ), which is defined in Sec.III [Eq.( 12)] and includes contributions only from paths with matching reactant and product energies, and which therefore resemble instanton-like configurations.The GR-QTST method for this one-dimensional system was evaluated using the same Monte Carlo importance sampling approach employed in Paper I. The results are, at least for low-dimensional systems, approximately independent of τ and so the exact specification of the value used to evaluate the GR-QTST rate is relatively unimportant.The ansatz of Paper I used a value of τ which maximizes the constrained action, φ c (τ ). 46However, in this paper, we choose to evaluate the rate at the value of τ = τ * which maximizes φ u (τ ).This is simpler to compute, especially when using the PIMD method described in Sec.III and the slight change of ansatz does not affect any of our conclusions from Paper I.Although GR-QTST is not rigorously separable for a system with two transition states, it is at least approximately so, as the rate is approximately independent of τ . 63For this reason, it is seen that GR-QTST accurately predicts the rates for this simple one-dimensional model consisting of two transition states.This is true whether the method is applied directly to both crossings simultaneously (GR-QTST) or whether the rate is obtained as the sum of the two individual rates (GR-QTST separated), as can be seen in Table I.In fact we showed in Paper I that GR-QTST reproduces the exact rate for a one-dimensional system of two crossed linear potentials at any temperature. 46Because in this special case, φ c is rigorously independent of τ , it would therefore also be exact for a system of many crossed linear potentials.
Note that there is a danger that for high-dimensional complex systems exhibiting multiple transition states, GR-QTST in its present form could suffer from similar problems to Wolynes theory, because it was explained in Paper I that φ c (τ ) becomes increasingly curved as the number of degrees of freedom is increased.It will therefore be important to address this size-consistency problem when applying it to large systems at low temperatures.However, one expects the problem to be less severe than for Wolynes theory because, unlike Wolynes theory, GR-QTST rigorously tends to the correct classical limit as α → 0, which is a general property of the method for all systems, irrespective of the number of degrees of freedom and number of transition states.This is because the GR-QTST method contains a constraint functional which forces the ring polymer (which is collapsed in this limit) to lie on diabatic crossing seam, thereby sampling only configurations at the classical transition state.

III. METHODS
Going beyond the simple system discussed in the previous section, it is important to discuss how GR-QTST rate can be computed for electron-transfer reactions in complex models and molecular systems, where the Monte Carlo importance sampling scheme described in Paper I falls short.Therefore, in this section we introduce a PIMD implementation of GR-QTST.In Sec.IV, we will then apply this method to an anharmonic and non-separable two-dimensional system and demonstrate that GR-QTST alleviates the problems identified with Wolynes and instanton theory.
In Paper I, we presented the formal theory of the GR-QTST method in terms of continuous path integrals which travel for imaginary-time τ on the reactant state, |0 , and for imaginary-time β −τ on the product state, |1 .Here we present only the working equations.
As in Wolynes theory, 9 the path integral is replaced by an ensemble average of discretized imaginary-time paths.These paths can be represented by a ring polymer of N ≡ N 0 + N 1 beads, in which N 0 − 1 beads feel the potential of the reactant state and N 1 − 1 feel that of the product state.The remaining two beads are called the "hopping beads" and feel a potential which is the average of the two surfaces.
We introduce λ as a dimensionless order parameter which determines how the beads are distributed between the two diabatic states and is defined by where because N 0 and N are integers, only discrete values are allowed.There are two special cases; when λ = 0, the ring polymer is completely on the reactant electronic state, and when λ = 1 it is completely on the product electronic state.
The unconstrained ensemble can be sampled using thermostatted PIMD 64 based on the following Hamiltonian: U (1) where T is the temperature, and x = {x (1) , . . ., x (N ) } are the positions of the beads with conjugate momenta p.The cyclic index i runs over each bead such that x (0) ≡ x (N ) , and the index j runs over each of the D nuclear degrees of freedom of the system.The reactant partition function is defined by an ensemble of ring polymers First we briefly discuss the computation of Wolynes rate theory, 9 given by: where F u (λ * ) is the maximum unconstrained free energy with respect to λ and thus also defines τ * through Eq. ( 6).It can be calculated using adiabatic switching and thermodynamic integration 66 where and . . .

(λ)
u denotes the canonical ensemble average with the Hamiltonian in Eq. ( 7), which can be obtained from PIMD simulations.
The GR-QTST rate formula, as mentioned in previous sections, is defined in terms of a constrained effective action, φ c (τ ), and the reactant partition function, Z 0 [Eq.( 8)].The constrained action is defined by and the definition of the function σ λ (x), which constrains the reactant and product energies to match, can be found in Appendix B. Direct calculation of the constrained action and the partition function individually is difficult using PIMD, but, as is typical in rate theory, only the constrained free energy, F c (λ), relative to the reactant is actually required.Although not calculated in this way, this free energy is defined by e −βFc(λ) ≡ e −φc(τ )/ /Z 0 .Thus the GR-QTST method approximates the golden-rule rate with 46 According to the ansatz of Paper I, λ * would be chosen as the value of the order parameter at which the constrained free energy is maximum.However, in this work, we use the simpler ansatz discussed in the previous section, in which λ * is the same as that used in Wolynes theory, i.e. the maximum of the unconstrained free energy.
In order to obtain the constrained free energy, we first express the delta function as the limit of a function of a new variable K, The introduction of this function makes it feasible to compute the GR-QTST rate with a PIMD free-energy simulation using the biased Hamiltonian One can see that this Hamiltonian includes an "umbrella-like" potential of the function σ λ (x).In contrast to the usual umbrella-sampling approaches, 66,67 here the centre of our umbrella is fixed but the strength can be varied.Note that we have chosen this definition such that when K = 0, the unconstrained ensemble will be obtained.However, as K tends to infinity, it enforces the GR-QTST energy constraint exactly. 68When performing the PIMD simulations with Hamilton's equations of motion based on H , an extra term appears in the first derivative due to the umbrella: An explicit formula for the derivative of σ λ is given in Appendix B.
Given that the unconstrained free energy F u (λ) has been calculated as described above [Eq.( 10)], we propose the following scheme for calculating the constrained free energy, 68 where and . . .
denotes the canonical ensemble average with the biased Hamiltonian in Eq. ( 15).Here we are using thermodynamic integration (TI) to measure the free-energy change upon slowing turning on the constraint and we will call this method for evaluating constrained free energies "δ-TI".
To avoid integration to infinity in Eq. ( 17) and to greatly improve numerical stability, the integral over K is performed using the coordinate transform where ξ ranges from 0 to 1 and K 0 is a scaling parameter, which we chose to be 0.1 in the following.The integral can now be performed using Tests showed that the results barely depend on K 0 values over two orders of magnitude, . We note that there are also other coordinate transformations that one could employ to perform the integral, but it was found that the above definition yields results reliably and efficiently for the systems studied here.
With the new ansatz, only one δ-TI needs to be performed to obtain the GR-QTST rate, which is for λ = λ * .This is simpler to compute compared with the GR-QTST ansatz of Paper I, for which several δ-TI calculations at different values of λ would be required.Our results will show that F c is approximately independent of λ for the systems investigated in this work, and thus the results are not significantly affected by this change.
Note that one could conceive of applying other free-energy sampling schemes to obtain the constrained free energy, which may have computational advantages.Possible alternative approaches which could be employed in future work include thermodynamic integration of the potential of mean force obtained from a constrained PIMD, metadynamics or adiabatic free-energy dynamics, 66 each using σ λ as the reaction coordinate.
The classical rate in the golden-rule limit is defined by 50 where Here F cl is the classical free energy of the crossing seam relative to the reactant well.It can be defined in a similar way to the ring-polymer versions given above with N = 1 and λ = 0 and can also be calculated with a δ-TI approach similar to Eq. ( 20).Finally, we shall also employ the semiclassical instanton method as described in Ref. 19 and the steepestdescent approximation (harmonic limit) to the classical rate, which can be calculated with the procedure described in Ref. 50.

IV. RESULTS
In the model system described in Sec.II, we were able to avoid the problems associated with Wolynes theory by computing the rate separately for the two product channels.
However, realistic systems are complex and multidimensional, and typically more than one transition state may exist for each product channel.Therefore performing Wolynes theory for each transition state separately will not be possible.In this section, we introduce a 2D model system in which there are two transition states, but only one product state.We show that Wolynes theory breaks down in the same way as shown in the 1D system but that nonetheless the GR-QTST method continues to give correct order-of-magnitude estimates for the rate constant when compared with numerically exact results.
Using the coordinates x = (x, y), the potentials of our 2D model are defined as where V 1 is chosen to be a separable system so as to simplify the computation of the exact rate.
However, the overall system is coupled and thus provides a nontrivial test for the methods.
The following set of parameters defines "system I": h 0 = 4.8, a 0 = 0.4, h 1 = 32, a 1 = 1, h 2 = 3.2, h 6 = 0.32, d = 0.75, p 1 = 28.8,p 2 = 1.15, p 3 = 0.7 (using reduced units with = 1).This system was designed to ensure that the gradient of V 0 is very different at the two transition states (they differ by a factor of 6.5), which means that the optimal τ values for the two crossings will be far apart.In Sec.IV D, we will tune some of the these parameters to create "system II".
A plot of the PES of system I is shown in Fig. 4(a).The potential mimics complex reactions in solution, where more than one TS exists due to the fluxional anharmonic solvent environment.We will study the electron-transfer rate for a particle with mass 35 at low temperature β = 3.Although these parameters were chosen to cause problems for Wolynes theory, we note that the breakdown of Wolynes theory occurs more generally, as was explained in Sec.II.To demonstrate this numerically, we also studied this system at higher temperature β = 1 where the reaction behaves classically.

A. Computational details
We computed rates in the golden-rule limit for system I using classical golden-rule transition-state theory, 5, 50 Wolynes theory, 9 semiclassical instanton theory, 18,19 and the newly developed GR-QTST method. 46Instanton theory was applied using the ring-polymer approach based on the Lagrangian formalism described in Ref. 19 with 64 beads on each PES.
In this method, the instanton configurations are obtained by a first-order saddle-point optimization in the combined space of x (bead positions) and τ variables.For comparison, numerically exact quantum results were calculated using a discrete variable representation (DVR). 69,70A grid of 200 × 50 DVR points was used, and we tested that a larger grid of 250 × 80 DVR points gives the same results.V 1 is separable in the x and y degrees of freedom, hence we used the same one-dimensional DVR as for V 0 in the y degree of freedom, and analytical scattering wave functions for the x degree of freedom. 71 performed Wolynes rate calculations with PIMD simulations using the Andersen thermostat and N = 128 ring-polymer beads to represent the imaginary-time paths.At β = 3, each PIMD simulation ran for 108,000 steps with a time-step of 0.03, with the first 8,000 discarded for equilibriation.For 56 ≤ N 0 ≤ 72, the simulations ran for 208,000 steps and employed replica exchange with 36 replicas (with the "hottest" replica at β = 0.105) to ensure ergodic sampling of the two crossings (see Appendix C for details).At β = 1, the simulation setup was almost the same as at β = 3, with the only difference in the replica-exchange setup.Here, replica-exchange sampling was performed only for 68 ≤ N 0 ≤ 96, with 16 replicas (the "hottest" replica at β = 0.11) and 308,000 simulation steps.The computed free-energy integrands are shown in the supporting information (SI). 72e GR-QTST simulations were also performed using N = 128.The simulation setup was mostly the same as the Wolynes simulations (108,000 steps), with two differences.First, a shorter time-step of 0.005 was used for simulations with K ≥ 1, because the constraining forces increase with the value of K.This does mean that the simulations become less efficient at large K, but nonetheless it was not a problem to converge the results for the systems investigated, such that we could obtain reasonably small error bars.Second, replica-exchange sampling was applied for all simulations, whatever the value of N 0 .The thermodynamic integration used at least 12 different K values in the range [0.0005, 16].The integral was computed with Simpson's rule, although we noted that using the trapezoidal rule instead gives the values within or comparable to the error bars of the results.
The classical rates were calculated with molecular dynamics simulations which ran for 205,000 steps with the first 5,000 discarded for equilibriation, at 17 different K values in the range [0.0001,16].The time-step and replica-exchange setup were the same as in the GR-QTST simulations.
The computational setup for system II was almost the same as for system I for all the methods, with the only difference in the replica-exchange setup.The Wolynes calculations on system II used replica exchange with 16 replicas (with the "hottest" replica at β = 0.15) for 42 ≤ N 0 ≤ 60, and 8 replicas for 60 ≤ N 0 ≤ 70, which ran for 208,000 steps.The GR-QTST calculation used 16 replicas with the "hottest" replica at β = 0.15, which ran for 108,000 steps.The free-energy integrands are shown in the SI.

B. Benchmarking the rates calculated using various methods
Here we discuss only the results obtained for system I, which are summarized in Table II, and leave the discussion of system II to Sec.IV D. There are two transition states (which we call "top" and "bottom") each with its own instanton, as shown in Fig. 4(a).The two transition states have a different character such that (at β = 3) the top instanton corresponds to λ = 0.29 whereas the bottom instanton corresponds to λ = 0.72.Instanton theory predicts that at this temperature, the reaction rate through either instanton is roughly TABLE II.Comparison of the golden-rule rates, k/∆ 2 , (with powers of 10 in parentheses) in reduced units for systems I and II, calculated with various methods.The classical rate is computed in two different ways: once in full using molecular dynamics (MD) and once using the steepestdescent (SD) approximation, which is equivalent to a local harmonic analysis around the TS.
The standard error in each simulation is estimated in the standard way taking into account the autocorrelation time. 67

System I
System II equal, as shown in Fig. 4(b).This occurs because although the bottom transition state is higher in energy, it has a thinner barrier which is more conducive to tunneling.
We find that the steepest-descent approximation does not lead to significant errors for this system.This is the case both for the harmonic classical rate, which is fairly close to the full classical result, as well as for the instanton, which is close to the exact quantum result.However, at β = 3, both classical methods underestimate the rate by more than an order of magnitude because nuclear quantum effects such as tunneling are not accounted for.
Instanton theory is thus seen to describe this tunneling effect very accurately, as expected from previous work, 20 and shows good agreement with the exact results, with only a 5% error mostly emanating from the harmonic approximation to the reactant partition function.
Wolynes theory however, overestimates the rate by two orders of magnitude, and actually log k has a larger error (albeit in the opposite direction) than the classical rate for this system.This example shows that the break-down of Wolynes theory identified in Sec.II is a general phenomenon, which may exist in many realistic reactions.
The GR-QTST rates are calculated using the δ-TI procedure described in Sec.III, and we show an example of the constrained free energy obtained and the thermodynamic integration over ξ in Fig. 5 (similar plots are given in the SI for the other systems).One can see that under the GR-QTST energy constraint, the free energy becomes almost independent of λ, implying that the GR-QTST rate will not be strongly affected by our choice of ansatz.The free-energy integrands show that the largest effect of the constraint comes from ξ < 0.5, corresponding to small K values (K < 0.1).The GR-QTST rate agrees well with the exact result, with only ∼15% error at β = 3 (which is comparable to the statistical uncertainty of the result).This confirms that the energy constraint cures the problem of Wolynes theory in these systems.At a higher temperature (β = 1), quantum tunneling becomes less important and the classical rates do not differ much from the exact quantum rate.Instanton theory also performs well as it is known that in this limit it tends to the harmonic classical rate, 18 which is a good approximation in this case.GR-QTST also rigorously tends to the correct classical limit, 46 predicting a rate within statistical error of the exact result.The Wolynes rate, on the other hand, does not tend to the correct limit, and it suffers from the same problems for this system as it did in the β = 3 case, overestimating the rate by a factor of 2.5 here.
These finding indicate that Wolynes theory could drastically overestimate the rate, seem-ingly predicting huge tunneling effects in systems that have multiple transition states, even in the high-temperature limit where a classical approximation would be valid.However, by including an energy constraint in the path-integral sampling, GR-QTST is able to give reliable predictions for the rate, both at high and low temperatures, with a similar accuracy to the rigorously-derived instanton theory.

C. Reaction mechanism
Obtaining an accurate prediction of the rate is only one objective for simulations of electron-transfer reactions, and characterizing the correct reaction mechanism is arguably even more important.Here we compare the mechanistic information obtained from Wolynes theory and from GR-QTST, based on the ensemble of ring polymers sampled by these methods.These should correspond to a transition-state ensemble and thus the relative populations at two transition states should quantify the "mechanistic branching fraction", which we will define as probability for a reactive particle to use one of the two available pathways.
We use the two instanton configurations and the values of two separate semiclassical instanton rate calculations to give a benchmark, i.e. k (s) SC /k SC , where here s labels either the top or bottom TS.
The unconstrained ensemble described in Sec.III with λ = λ * is used for a mechanistic analysis of Wolynes theory.On the other hand, the relevant ring-polymer ensemble for GR-QTST is obtained by performing constrained PIMD simulations (also with λ = λ * ) with the constraint σ = 0, using the RATTLE algorithm, 73 and re-weighting to recover the soft constraint limit. 66,67Assuming that the ensemble makes a clear split into two separate regions of configuration space, the probabilities corresponding to the mechanistic branching fraction can be extracted.
The density plots of the ensembles sampled from constrained and unconstrained simulations are shown in Fig. 4(a).One sees immediately that the transition-state ensemble according to Wolynes theory does not include the top instanton path and will therefore give an unphysical interpretation of the mechanism of this reaction. 74On the other hand, in GR-QTST, the energy constraint ensures that the samples remain close to the crossing seam and thus closer to the instantons in the system, predicting the correct reaction pathways.The branching fraction for the mechanism through the bottom TS is given by instanton theory as 32%.A similar result of 39% is found by GR-QTST and 44% by Wolynes theory.Therefore, despite the fact that Wolynes theory overpredicts the rate by two orders of magnitude, the prediction for the branching fraction is seemingly not so bad.As we will explain shortly, this is a coincidence due to the fact that the two pathways have a roughly equal contribution.
In the higher-temperature case (β = 1) where the Wolynes theory predicts a more reasonable reaction rate (within an order-of-magnitude of the exact result), the predicted reaction mechanism is, however, qualitatively incorrect.As shown in Fig. 4(b), instanton theory predicts that at this temperature transitions through the top transition state is the major contribution to the rate, whilst only 3% of the reactive particles will pass through the higher-energy bottom transition state.As shown in Fig. 6, Wolynes theory yields an incorrect physical picture, in which both TSs have a significant contribution to the rate at β = 1, with the bottom TS contributing 36% to the reaction.In contrast, GR-QTST identifies that only the top TS is important, and that the bottom TS has only a small contribution of 3%, giving the correct physical picture, in agreement with the rigorously-derived instanton approach.
In fact, whenever Wolynes theory breaks down for systems with two TSs, it will typically sample both TSs as if they are equally important pathways and give a mechanistic branching fraction of roughly 50%, regardless of the true dynamics.We can justify this statement based on Fig. 1 and the discussion in Sec.II which suggest that τ * will typically be found approximately at the location where the two constituent φ (s) u (τ ) terms are equal to each other and thus give roughly equal contributions to the rate.

D. System with strong anharmonic fluctuations
The GR-QTST method was developed to go beyond the semiclassical approximation of instanton theory.The accuracy of the semiclassical approximation depends on the anharmonicity of the fluctuations around the instanton path.In this section, we design a strongly anharmonic system, for which instanton theory leads to a significant error, and test the accuracy of the GR-QTST against exact benchmark results.
We modified some of the parameters of system I to create system II.In particular, we defined d = 0.9 and then chose new values of three parameters h 2 = 4.848, p 2 = 0.7145, p 3 = 0.59 to give an almost zero-frequency normal mode of the instanton at the temper- ature β = 3. 75 The existence of this very low frequency causes higher-order fluctuations to dominate the path integral, which are neglected by the semiclassical approximation.A contour plot of system II together with the instanton is shown in Fig. 7.We note that, with system II in particular, it would have been very inefficient to use the importance sampling approach based on normal modes as described in our previous work. 46However, the new PIMD implementation allows such anharmonic systems to be simulated with a reasonable computational effort.
The golden-rule rate constants calculated with various methods are given in the final column of Table II.The classical rates are more than two orders of magnitude lower than the exact benchmark, suggesting that this reaction is in the deep tunneling regime where nuclear quantum effects are prominent.Note that there is close agreement between the classical harmonic approximation and the full classical rate calculated with molecular dynamics sampling.This is because the classical transition state does not have frequencies close to zero and thus anharmonic effects are not so important here.Instanton theory, on the other hand, overpredicts the rate by almost an order of magnitude.This is because the neglect of higher-order fluctuations causes it to overestimate the true impact of fluctuations around the path.The Wolynes rate is ∼70% too high, because it suffers from the same problems discussed in the previous section, albeit to a lesser extent for this set of parameters.
The GR-QTST method makes neither the saddle-point approximation nor the semiclassical approximation, and because of this, is seen to predict a rate much closer to the exact results.

V. CONCLUSIONS
In this work, we have presented an approach for computing the GR-QTST rate from PIMD simulations, allowing the method to be applied efficiently to multidimensional anharmonic systems.We have applied this method to a model system with two transition states and demonstrated that GR-QTST gives accurate predictions for the golden-rule rate constant, whereas Wolynes theory can drastically overestimate the rate.This error of Wolynes theory can occur even in the high-temperature regime, such that there is a danger of predicting strong quantum tunneling effects where there are none.We explained that the cause of this problem is the saddle-point approximation employed by Wolynes theory.The energy constraint in the GR-QTST method ensures that paths are sampled near to the instantons of this system, whereas Wolynes theory does not, and it cannot therefore be used to extract mechanistic information.The GR-QTST method also predicts accurate rates for a system where semiclassical instanton theory breaks down due to its neglect of anharmonicity in the fluctuations around the instanton.Together, these promising results indicate that, at least for certain types of systems, the GR-QTST method is the most accurate imaginary-time path-integral approach for calculating electron-transfer rates in the golden-rule limit.
The GR-QTST method does however still require further development and testing.In Paper I, 46 we explained that the method, as it is currently defined, is not rigorously size consistent in general, although it nonetheless was found to give accurate results for a multidimensional spin-boson model and can be shown to rigorously give the correct result in the classical limit for any system.In the future, we will test whether GR-QTST can maintain this good performance when applied to atomistic simulations of solvated reactions.surfaces and x 0 is the location of both transition states, which have the same energy by construction.
We perform calculations in reduced-dimensional units.The degree of the quantumness is defined by α = β ω, the dimensionless barrier height by Φ = βmω 2 x 2 0 /2 and the reduced gradient by η s = − κs mω 2 x 0 .At first, we study the system in the quantum regime where tunneling is significant.We have chosen α = 2.5, Φ = 45, η A = 0.5 and η B = 2 to obtain the results in Table I. Results are presented as the dimensionless ratio k/k cl , where The path integral representation of Eq. ( 3) is Gaussian and can be evaluated analytically such that the exact effective action for the system of Eq. (A1) can be written in dimensionless units as follows: e −φ (s) u (τ )/ = Z ‡ (τ ) e − S(s) (τ )/ , (A2) where the classical action S(s) (τ ) is and the effective transition-state partition function is These formulae can be used in Eq. ( 2) with complex argument, τ − it, to obtain the exact result by numerical integration over time.Alternatively, they can be used with t = 0 to evaluate the result of Wolynes theory, Eq. ( 4), by numerically solving for the optimal τ * .
The reactant partition function in reduced dimensions is simply Z 0 = [2 sinh( 1 2 α)] −1 and is not affected by the semiclassical approximation.
The semiclassical approximation is very accurate for this system, and the following formula was used to compute the instanton rate:  This expression is evaluated at τ = τs which satisfies the equality d S(s) dτ = 0. We can study the expression in certain limiting cases similarly to what was done with the instanton rate for Ĥ0 ] is the reactant partition function and the unconstrained effective action, φ u (τ − it), is defined by the correlation function e −φu(τ −it)/ = Tr n [e − Ĥ0 (τ −it)/ e − Ĥ1 (β −τ +it)/ ].

2 FIG. 4 .
FIG. 4. (a) Contour plot of the two potential-energy surfaces of system I.The blue-green contours show V 0 and the red-yellow contours show V 1 .Only the lower of the two potential-energy surfaces is shown at each coordinate and the dashed line shows the crossing seam between the two surfaces.The '×'s mark the two classical transition states on the crossing seam and the two instantons for β = 3 are also shown, with blue beads on V 0 , red on on V 1 , whereas the violet bead indicates the hopping point.The filled contours represent density plots of the bead positions at β = 3 for unconstrained/Wolynes (grey) and constrained/GR-QTST (green) obtained from simulations with λ = λ * .(b) A representation of the "mechanistic branching fraction," defined in terms of contributions to the total rate through the two TSs, here according to instanton theory at β = 3 and β = 1.

75 FIG. 5 .
FIG. 5. Free-energy calculations.(a) Comparison of the λ dependence of the unconstrained free energy (used for Wolynes theory) and the constrained free energy (used for GR-QTST) for system I at β = 3.The dashed arrows illustrate how the GR-QTST free energies are obtained (Eq.(17)).(b) The constrained free-energy integrand at different λ values obtained from biased PIMD simulations.The colours correspond to the colours of the dashed arrows in (a).The statistical error bars in this case are smaller than the symbol size.

FIG. 6 .
FIG. 6.Comparison of the distribution of y + variables, defined as the y-component of the average hopping point (see Appendix B), sampled by the unconstrained ensemble of Wolynes theory and the constrained ensemble of GR-QTST, both at λ = λ * with β = 1.

2 FIG. 7 .
FIG. 7. Contour plot (as in Fig. 4) of system II, which has only one TS.Here only one instanton exists at β = 3. Due to the strong anharmonicity the instanton has clearly moved far from the TS due to corner-cutting effects.The filled contours show density plots of the bead positions at β = 3 for unconstrained/Wolynes (grey) and constrained/GR-QTST (green) ensembles obtained from simulations with λ = λ * .