- Email: [email protected]
A dynamical test of the centroid formulation of quantum transition state theory for electron transfer reactions
Volume 206, number 1,2,3,4
CHEMICALPHYSICSLETTERS
30 April I993
A dynamical test of the centroid formulation of quantum transition state theory for electron transfer reactions C.H. Mak Departmentof Chemistry,Universityof SouthernCalifornia, Los Angeles,CA 90089-0482,USA
and John N. Gehlen Departmentof Chemistry,Universityof Caltjornia,Berkeley, CA94720, USA Received 8 December 1992;in final fonn 11 February 1993
Rigorous quantum dynamics simulations have been performed to test the path-integralcentroid formulation proposed by Voth, Chandler and Miller for quanta1rate constants in a simple class of models for electrontransferreactionsin the tight-binding limit. Fully dynamical quantum Monte Carlo simulations show that inside and slightlyaway from the diabatic limit the electron flux is invariant with the electronic centroid density near the dividing surface instead of being sharply peaked there. Despite this, we show analyticallythat the rate in this regime is largelygoverned by the ccntroid density on the dividing surface. But on approaching the adiabatic limit with largerelectronic coupling,complex large-amplitudecoherent oscillations develop in the quantum flux, suggestingthat in addition to the imaginary-time centroid density, the dynamical factors also play an important role in determining the full quantal electron transfer rate.
1. Introduction In the theory of classical activated processes, the reaction rate constant can be rigorously factored into two factors [l-4]: (a) the transition state theory estimate for the rate constant based entirely on an equilibrium theory, and (b) a dynamical correction factor (sometimes called the transmission coeficient) which accounts for dynamical recrossing events. This rigorous factorization forms the basis for the exact calculation of reaction rates and at the same time provides a computationally feasible procedure for treating classical rare events [ l-51. Since the rigorous statistical mechanical derivation of the exact classicalrate expression and the development of the corresponding computational methods [l-6], the search for a similar correlation function formulation and similar computational methods for treating rare events in quantum systems has produced little success. But recently, Voth, Chandler and Miller [ 7 ] (hereafter referred to as
130
VCM) presented a correlation function analysis for the quantum rate constant and based on a path integral formulation of the rate expression, they suggested an interesting factorization of the quantum rate constant similar to the classical theory. Their formulation was motivated by Gillan’s earlier paper [ 8 1. In that work, it was observed that a computational method based entirely upon imaginary-time path integrals reproduces both the high-temperature classical Arrhenius factor and the low-temperature instanton expression for the rate constant in a double-wellpotential linearly coupled to a classicalharmonic bath, although a different velocity factor was required in each limit. Adopting Gillan’s idea that the quantum rate constant is related to the equilibrium statistics of the centroid of the imaginary-time quantum paths, VCM formulated rigorous expressions for the quantum transition state estimate for the tunneling rate, as well as the dynamical frequency factor. At about the same time as Gillan, a related imaginary-time path integral approach was
0009-2614/93/S 06.00 Q 1993 Elsevier Science Publishers B.V. AUrights reserved.
CHEMICAL PHYSICS LETTERS
Volume 206, number 1,2,3,4
proposed by Wolynes for electron transfer rates in the Golden Rule limit [9], and this approach was later recast in terms of the centroid formulation by Chandler [ lo] and Bader, Kuharski and Chandler
[illIn this Letter, we present results from what we believe is the first quantum Monte Carlo simulations of the electron transfer rate constant in the tightbinding limit using real-time path integrals and compare them to VCM’s centroid formulation of the quantum rate constant. These simulations are exact numerical reproductions of the fully quantum-mechanical Feynman integral representation for the dynamics of the electron plus an infinite number of bath modes in the deep-tunnelingregime. Our results show that the dominant contribution to the quantum rate constant in the Golden Rule limit is indeed captured by the equilibrium statistics of the centroid of the imaginary-time quantum paths, but for a rather different reason from that originallysuggestedby VCM. We will examine these findings in the Golden Rule limit in light of an analytic continuation theory. For systems outside the Golden Rule regime, we show that the simulated quantum flux exhibits complex quantum behaviors which are not well understood within imaginary-time theories.
2. The centroid formulationof quantumtransition state theory We briefly review VCM’scentroid formulation for the quantum reaction rate to establish the notations used in the rest of this Letter. The starting point is the quantum rate expression derived from correlation functions [ 121: k A_B=k(t)=
30 April1993
real- and imaginary-time propagators in eq. ( 1) are replaced by functional integrals over all paths. Fig. 1 shows a diagrammatic representation of the integrals involved in the calculation of k,,, [ 121.In fig. 1, the hypersurface 8 divides space into A and B states. The solid line represents the action for (41 exp( -BH) 14'), and the dashed and dotted lines the actions for (q’ 1exp(i&/h) (4”) and (4” /exp( i%/A) 1q). The integral runs over all positions q on the A side, all positions q’and q’ on the B side, and finally over all possible paths connecting these three vertices. By a coordinate transformation, one can project out the centroid of the thermal path qo=(@)-’ X16 dtq(t), where q(0) corresponds to q in fig. 1 and q(@i) to q’. VCM factor eq. ( 1) in such a way to focus attention on the role of the thermal path centroid q,,
where p(q)=Q-'l~q(t)d(qo-q)exp{-S[q(t)l} is the equilibrium density for finding the centroid at position q and S is the action for the thermal path, v(q; T) is a dynamical frequency factor which encompassesall the remaining path integralsin eq. ( I ) , and I is defined by the last equality. Apart from a multiplicative constant, the product p(q)v(q; F) is equal to the quantum flux across the dividing surface with the centroid of the thermal path constrained at position q. VCM argue based on analytical considerations that the dynamical frequency factor y( q; f)
-? @ix* Q-’
xImTr[hAexp(-/?H) x h, exp(itH/fi) he exp( -ilH/fi)]
,
(1)
where hA and hg= 1-hA are the population operators of the reactant state and the product state, respectively, H is the Hamiltonian, Q is the canonical partitionfunction,x,=Q-‘Tr[h,exp( -pH)] isthe equilibrium mole fraction of state A and tis the plateau time [ 11. Using Feynman’s path integrals, the
Fig. 1. Diagram representing the path integral in eq. ( 1). The solid line represents an imaginary-time path; the dashed and dotted lines represent real-time paths. The surface @ divides space into A and B states. The integral runs over all positions 4 on the A side, all positions q’and q”on the B side and over all paths.
131
Volume 206, number I ,2,3,4
CHEMICAL.PHYSICSLETTERS
should be sharply peaked at @ such that the dominant contribution to the total flux will come from q near the dividing surface, and the factorization iu the second equality in eq. (2) has been made to emphasize this.
3. Hamiltonian model for electron transfer reactions In the following, we study the nature of the centroid density p(q) and the dynamical frequency factor u(q; f) of electron transfer reactions in the tightbinding limit within a class of simple Hamiltonian models, H=H,
tHB(X) tp[ B(X)t Jo] )
1 dtC(t), --m
W
where C(t)=Q-‘Tr{exp[ xexp( -itH,/h)}
-(fi-it)H,/h] ,
(4b)
and HA and Hu are the diabatic (K=O) Hamiltonians of states A and B. Wolynes [ 91 suggestedthat the integral in eq. (4a) could be more efficiently computed by distorting the integration contour to go from t* -co to r? t co such that the contour will pass through a stationary point t* in the complex-t plane. For the Hamiltonian in eq. (3) with &,=O, a sta132
tionary point does exist at t* = - f ifi, and as has been pointed out in refs. [7,10,11], the value of the integraud C( t* ) is precisely the centroid density p ($ ) on the dividing surface. Therefore, the centroid formulation in the Golden Rule limit appears to correspond exactly to the saddle point approach of Wolynes. Notice, however, that the reaction coordinate considered in these previous studies aud the present one is a discrete variable in the tight-binding limit, whereas the VCM centroid formulation was originally proposed for a continuous reaction coordinate.
4. Quantum dynamics simulations of electron transfer rate
(3)
where Hm= -K&; and & is the Pauli spin operator, HB is a bath Hamiltonian, ,L=,ut$ is the dipole operator, b= Ci cix, is the solvent polarization with xi being the coordinate of the ith bath mode, and &#O for asymmetric electron transfer. In the following, we will treat the symmetriccase &o=0 only. We also assume HB is harmonic, in which case the influence of the bath on the two-level system is captured entirely by the bath spectral density J(o) [ 13,141, which is related to the bath correlation function ( B( 0) J(t) ) [ lo]. Notice that within this tight-binding model, the position of the electron is described by a discrete spin variable 0 - the electron is in the A state if 6= + 1 and in the B state if cr= - 1. In the limit K-+0, the rate can be calculated using the Golden Rule of perturbation theory [ 9,101: xAk~~B=(K/fi)z
30 April 1993
To understand the factorization made in the centroid formulation when it is applied to the tightbinding electron transfer model, we study in great detail the behavior of the centroid density p(q) of the electronic’coordinate and the dynamical frequency factor v(q; t) as a function of q with rigorous quantum Monte Carlo methods. We adopt the aqueous Fe2+/Fe3+ self-exchange reaction as a model, using the Hamiltonian in eq. (3) with the spectral density J(o) derived from the molecular dynamics simulationsof Kuharski et al. [ 15] at room temperature. With the classical isomorphism [ 161, the quantum system is mapped onto a 1D king model which has nonlocal influence functional bonds [ 16181. The calculations of p(q) involve purely imaginary-time path integrals and were carried out using standard umbrella sampling methods. The real-time quantum Monte Carlo simulations of v(q; t ) were performed using a filtering scheme proposed by Mak [ 191 in which the natural weight of each path q(t) in state space is filtered by a damping function D[q(t)]=l+{(1)[q(t)]+<(2)[q(t)], where {(I) [q(t) ] and c(*) [q(t) ] are the sum of the weights of all first-order and second-order fluctuations from the path q(t) defined in ref. [ 191 #I. xl Since the centroid of the imaginary-time path is constrained, not all first- and second-order fluctuations are allowed. The results reported here were obtained using all first- and second-order fluctuations and are therefore subject to a small uncertainty in the centroid position. However, the error introduced is small as long as the number of discretizations in the imaginary-time path is large,andcan bc further reducedby increasingthis number.
CHEMICALPHYSICSLETTERS
Volume 206, number 1,2,3,4
-1
-a.5
0
0.5
1
0
Q
0.2
0.4
30 April 1993
-1
-0.5
Kt
0
0.5
1
9
Fig.2. MonteCarlo results for,W=0.4. (a) Centroid freeenergycomputedusing 100 discretizations for the thermal path. (b) Dynamical frequency factor v(q;t)for several values of q.From top to bottom, q= 0,O.1250.25 and 0.375. 32 discretizations were used for the thermal path and each of the real-time paths. v(q; t) reaches a plateau after a short initial transient. (c) The plateau value v(q; r) as a function of q.
We first examine the results for small electronic coupling jlK= 0.4, which is within the Golden Rule limit. Fig. 2a shows the centroid free energy BF( q) = - In p ( q). As expected, the centroid free energy has a maximum at the dividing surface q= @ = 0, so the centroid density p(q) would have a pronounced minimum at @*.The dynamical frequency factors v(q; t) obtained from quantum dynamics simulations are shown in fig. 2b for several values of q. After an initial transient, v (q; I) for all q approach their plateau values v(q; L). Fig. 2c shows that the plateau values v( 4; t) decay rapidly as a function of q away from the maximum at g=@, and as VCM suggest, v (q; f) is indeed sharply peaked at 8. According to eq. (2), the rate does not depend on either p( q) or v( q; f) individually, but on the product p(q)v(q; f), which is proportional to the quantum flux with the thermal path centroid constrained at q. VCM suggest that because the dynamical frequency factor is sharply peaked at @, the flux should be sharply peaked there as well. As VCM demonstrate, this behavior is certainly true for classical or nearly classical systems [ 71. These are systems with continuous potentials along the reaction coordinate and above the crossover temperature T,: defined by AcoJkT, < 214 where U, is the barrier curvature. But for the highly quanta1 electron transfer system studied herein, the numerical results derived from the data in figs. 2a and 2c and plotted in f%. 3 show a different behavior, They show that the flux is constant within statistical errors for a broad range of q. The sharp decrease of Y(q; f) away from #’ appears to be delicately balanced by the rapidly increasing
p(q) so that the product is invariant with q. Because of the rapidly decreasing magnitude for v( q; f) away from q* = 0, substantially longer sampling times are required to compute v(q; f) with the same relative error. The time required to obtain the last point in fig. 3, which corresponds to q= to.4375 was a factor of 100 longer than q=O, and to further reduce the error to the same level as q=O would require a factor of 2500 more computational effort. To extend our calculations beyond the present limit would require impractically long sampling times. The results in fig. 3 strongly suggest that in the Golden Rule limit, p( q)v( q; f) is invariant with q in the tight-binding limit and is not dominated by the flux at a single point q*, in contrast with VCM’s observation for a continuous reaction coordinate. However, VCM’s suggestion that the rate is ultimately governed by the value of p( q) at q= fl is still
-I
-0.5
0
0.5
I
9
Fig. 3. The quantum fluxp(q)u(q; T) normalized top(q)v(q;i) plotted versus the centroid position q for /X=0.4. Vertical bars indicate one standard deviation. Dashed line denotes Golden Rule limit.
133
Volume206,number1,2,3,4
CHEMlCAL PHYSICS LETTERS
correct in the tight-binding model, at least in the Golden Rule limit. This is because if p( q) v (q; fl is indeed constant throughout the whole range of g, then according to eq. ( 2 ) /c*_~ is simply p (8) multiplied by a constant f I where r= 4 v(8; f). In this case, a short real-time quantum Monte Carlo simulation used to compute v(q; F) at a single point q= 8 plus a straightforwardimaginary-timecalculationfor p( q) are sufficient to compute an exact quantum rate. Having demonstrated numerically that in the Golden Rule limit p(q)v(q; f) is constant, we will now give a mathematical proof for it. The Golden Rule expressionin eq. (4) can be derived easily from eq. ( 1) by retaining only those diagrams that have exactly two kinks - one in the imaginary-time path and the other in one of the real-time paths [ 121. The kink in the imaginary-time path may occur at an imaginary time t anywhere between 0 and m so we can classify contributions to the Golden Rule rate constant according to the location of the kink along the imaginary-time path, (sa) where k”‘(~)=(K/h) xTr{exp[ -
7 dtQ-’ -co (@-7)HAlfi]
xexp( -it&/fi)
exp(itWh)
exp( -.r&)}
.
(5b)
Analytically continuing the integrand into the complex-t plane and using Cauchy’s theorem, it is easily shown that kc’) (7) is invariant with z. Within the two-state model defined in eq. (3)) z is related to the centroid of the imaginary-time path q. by a simple relationship #2:qo= 1- 2~/Bjl.It follows that k(‘) (r) is identical to p( qo)v ( qo; f) apart from a multiplicative constant, and since kc2’(z) is z invariant, p(qo)v ( qo; r) is also invariant with qo. These arguments also establish a clear relationship between a1This is because the time 7 separates the imaginary-time path into two segments, one having all down spins and the other hav-
ing all up spins. The length ofthe first segment is c and the second is J%I- 7. Hence the centroid of the whole imaginary-time path is so=(P)-‘I-r+(Bh-r)l.
134
30April 1993
Wolynes’imaginary-time path integral approach [ 91 and the centroid formulation in refs. [7,10,11]. Having demonstrated that the flux is invariant with the position of the centroid in the Golden Rule limit, we now return to the quantum Monte Carlo results and examine the behavior of p (q) and v(q; f) outside the Golden Rule limit. Fig. 4a shows Monte Carlo results for Y(q; t ) as a function of t for a larger electronic coupling j?K=2, which is well outside the Golden Rule regime. These results and data from centroid density calculations are summarized in table 1. The centroid density p (q) = exp [ - BF( q) ] and the dynamical factor v(q; f) display behaviors similar to the data shown earlier in fig. 2 for /X=0.4. Within statistical error, the quantum flux p( q)v( q) for /X=2 is approximately invariant with the centroidpositionqforO
0
0 0
i Kt
2
0
L
2 Xl
Fig. 4. Monte Carlo results for the dynamical frequency factor: (a) j3K=2, q=O, 0.125,0.25 and 0.375 from top to bottom. 32 discretizations were used for the thermal path and each of the real-time paths; (b) /WA 5, q=O, 0.125 andO.25 from top tobottom. 48 discretizations were used for the thermal path and 24 for each of the real-time paths.
Volume 206,number1,2,3,4
CHEMICALPHYSICS LETTERS
30 April 1993
Table 1 Results of Monte Carlo simulations for flK= 2
ww
4
4%
0 0.125 0.25 0.375
1 OS967 0.1516 0.01323
T)
P(q)
P(f7MK i)/P(W(o; i)
Error (%)
MC passes
0 0.4546 1.8314 4.1516
1 0.9401 0.9464 0.8406
0.71 0.69 2.82 12.5
0.5x 106 1.5 2.5 24.0
exhibits the usual plateau after an initial transient. However, for 4 close to the dividing surface, v (q; t ) displays curious oscillations which appear to persist for rather long times. These oscillations are the signature of coherent quantum processes. Since K determines the intrinsic tunneling frequency of the charge, if K is much larger than the typical bath frequency, the charge cannot be fully solvated, and coherence can persist for an extended duration. But interestingly, the results in fig. 4b also suggestthat the coherence appears to be damped rapidly as the centroid moves away from the dividing surface towards the stable states. The appearance of coherence for large values of K creates serious problems for the direct application of the centroid formulation to simulation of quantum electron transfer rates in the crossover region between the diabatic and adiabatic limits. Since largeamplitude coherent oscillations in the flux may persist for an extended period of time, a long real-time path integral calculation may be necessary to reach the plateau time, and as a result, such quantum simulations may be impractical. At the same time, the appearance of coherent motion implies the prevalence of kink formation within the real-time paths. Therefore, it is not clear whether a method based solely on imaginary-time calculations alone can reproduce an electron transfer rate which is possibly highly influenced by coherent processes. Our results suggestthat the centroid breakup of the full electron transfer rate indicated by eq. (2) in terms of the electronic centroid is perhaps not optimal for systems approaching the adiabatic limit. Based on the physical situation, it may be argued that a similar breakup of the electron transfer rate in terms of the nuclear centroid is more appropriate [20], but the validity of this alternative formulation remains to be tested [2 11. We should point out that although our simulations suggestthat the centroid formulation for
the electron transfer rate using the electronic centroid may not be optimal for systems in between the diabatic and adiabatic limits, a recent calculation [ 221 of the electron transfer rate based on the electronic centroid free energy alone yields an activation energy that appears to be consistent with the conventional estimates in both the diabatic and adiabatic limits [24]. Finally, we also note that for strictly adiabatic systems, a computational scheme based on an entirely classical treatment of the solvent modes is sufficient for computing the electron transfer rate constant [ 231 #3.
Acknowledgement Throughout this work, the authors have benefited from many enlightening discussions with David Chandler, Masako Takasu, Reinhold Egger,Gregory Voth and Joel Bader. CHM wishes to thank David Chandler for introducing him to this area and for suggestingthe studies reported here. We also thank David Chandler and Greg Voth for helpful comments on the manuscript. This research was supported by the Camille and Henry Dreyfus Foundation, the USC FRIF Program and the Army Research Office. Computational resources from IBM are gratefully acknowledged. (iaFor a recent application to an electron transfer reaction in the strictly adiabatic limit, see ref. [ 241. For a related treatment of adiabatic proton transfer, see ref. [ 251.
References [ 1I D.Chandler, J. Chem. Phys. 68 ( 1978) 2959. [2 J J.T. Hynes,irx The thcoq of chemicalreactions,ed. M. Baer (CRC Press, Bcca Raton, 1985).
135
Volume 206, number 1,2,3,4
CHEMICALPHYSICSLETTERS
[ 31 B.J. Beme, in: Multiple time scales, eds. J.V. Brackbilland 9.1. Cohen (Academic Press, New York, 1985) . [ 41 R.O. Rosenberg, B.J. Beme and D. Chandler, Chem. Phys. Letters 75 (1980) 162. [ 5] C.H. Bennett,in: Algorithmsfor chemicalcomputation, ACS Symp.Ser. No. 46, ed. R.E. Cbristofferson (Am. Chem. Sot., Washington, 1977). [6] T. Yamamoto, J. Chem. Phys. 33 (1960) 281. [ 71 G.A. Voth, D. Chandler and W.H. Miller, J. Chem. Phys. 91 (1989) 7749. [8] M.J. Gillan, J. Phys. C 20 (1987) 3621. [ 9) P.G. Wolynes,J. Chem. Phys.,87 ( 1987) 6559. [lo] D. Chandler, in: Liquids, freezing and the glass transition, Les Houches 51, Part 1, eds.D. Levesque, J.P. Hansen and J. Zinn-Justin (Elsevier, Amsterdam, 1991). [ 111J.S. Bader, R.A. Kuharshi and D. Chandler, J. Chem. Phys. 93 ( 1990) 230. [ 121G.A. Voth, D. Chandler and W.H. Miller, J. Phys. Chem. 93 (1989) 7009. [ 131A.O.Caldeiraand A.J. Leggett,Phys. Rev. Letters 46 ( 1931) 211.
136
30 April 1993
[ 141A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fischer, A. Garg and W. Zwerger, Rev. Mod. Phys. 59 (1987) 1. [ 151R.A. Kuharski,J.S. Bader,C. Chandler,M. Sprik,ML. Klein and R.W. Impey, J. Chem. Phys. 89 (1988) 3248. [ 161D. Chandler and P.G. Wolynes,J. Chem. Phys. 74 (1981) 4078. [ 171M. Suzuki,Pros. Theoret. Phys. 56 (1976) 1454. [ 181D. Chandler, Introduction to modem statistical mechanics (Oxford Univ. Press, Oxford, 1987) [ 191C.H. Mak, Phys. Rev. Letters 68 (1992) 899. [ 201D. Chandler, private communication. [ 211C.H. Mak, in preparation. [ 22 1J.N. Gehlen, D. Chandler,H.J. Kim and J.T. Hynes, J. Phys. Chem. 96 ( 1992) 1748. [ 231R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 8 11 (1985)265. [ 241D.A. Zichi, G. Ciccotti, J.T. Hynes and M. Ferrario, J. Phys. Chem. 93 (1989) 6261. [ZS] D. Borgis, G. Tarjus and H. Azzouz, J. Phys. Chem. 96 (1992) 3188.
CHEMICALPHYSICSLETTERS
30 April I993
A dynamical test of the centroid formulation of quantum transition state theory for electron transfer reactions C.H. Mak Departmentof Chemistry,Universityof SouthernCalifornia, Los Angeles,CA 90089-0482,USA
and John N. Gehlen Departmentof Chemistry,Universityof Caltjornia,Berkeley, CA94720, USA Received 8 December 1992;in final fonn 11 February 1993
Rigorous quantum dynamics simulations have been performed to test the path-integralcentroid formulation proposed by Voth, Chandler and Miller for quanta1rate constants in a simple class of models for electrontransferreactionsin the tight-binding limit. Fully dynamical quantum Monte Carlo simulations show that inside and slightlyaway from the diabatic limit the electron flux is invariant with the electronic centroid density near the dividing surface instead of being sharply peaked there. Despite this, we show analyticallythat the rate in this regime is largelygoverned by the ccntroid density on the dividing surface. But on approaching the adiabatic limit with largerelectronic coupling,complex large-amplitudecoherent oscillations develop in the quantum flux, suggestingthat in addition to the imaginary-time centroid density, the dynamical factors also play an important role in determining the full quantal electron transfer rate.
1. Introduction In the theory of classical activated processes, the reaction rate constant can be rigorously factored into two factors [l-4]: (a) the transition state theory estimate for the rate constant based entirely on an equilibrium theory, and (b) a dynamical correction factor (sometimes called the transmission coeficient) which accounts for dynamical recrossing events. This rigorous factorization forms the basis for the exact calculation of reaction rates and at the same time provides a computationally feasible procedure for treating classical rare events [ l-51. Since the rigorous statistical mechanical derivation of the exact classicalrate expression and the development of the corresponding computational methods [l-6], the search for a similar correlation function formulation and similar computational methods for treating rare events in quantum systems has produced little success. But recently, Voth, Chandler and Miller [ 7 ] (hereafter referred to as
130
VCM) presented a correlation function analysis for the quantum rate constant and based on a path integral formulation of the rate expression, they suggested an interesting factorization of the quantum rate constant similar to the classical theory. Their formulation was motivated by Gillan’s earlier paper [ 8 1. In that work, it was observed that a computational method based entirely upon imaginary-time path integrals reproduces both the high-temperature classical Arrhenius factor and the low-temperature instanton expression for the rate constant in a double-wellpotential linearly coupled to a classicalharmonic bath, although a different velocity factor was required in each limit. Adopting Gillan’s idea that the quantum rate constant is related to the equilibrium statistics of the centroid of the imaginary-time quantum paths, VCM formulated rigorous expressions for the quantum transition state estimate for the tunneling rate, as well as the dynamical frequency factor. At about the same time as Gillan, a related imaginary-time path integral approach was
0009-2614/93/S 06.00 Q 1993 Elsevier Science Publishers B.V. AUrights reserved.
CHEMICAL PHYSICS LETTERS
Volume 206, number 1,2,3,4
proposed by Wolynes for electron transfer rates in the Golden Rule limit [9], and this approach was later recast in terms of the centroid formulation by Chandler [ lo] and Bader, Kuharski and Chandler
[illIn this Letter, we present results from what we believe is the first quantum Monte Carlo simulations of the electron transfer rate constant in the tightbinding limit using real-time path integrals and compare them to VCM’s centroid formulation of the quantum rate constant. These simulations are exact numerical reproductions of the fully quantum-mechanical Feynman integral representation for the dynamics of the electron plus an infinite number of bath modes in the deep-tunnelingregime. Our results show that the dominant contribution to the quantum rate constant in the Golden Rule limit is indeed captured by the equilibrium statistics of the centroid of the imaginary-time quantum paths, but for a rather different reason from that originallysuggestedby VCM. We will examine these findings in the Golden Rule limit in light of an analytic continuation theory. For systems outside the Golden Rule regime, we show that the simulated quantum flux exhibits complex quantum behaviors which are not well understood within imaginary-time theories.
2. The centroid formulationof quantumtransition state theory We briefly review VCM’scentroid formulation for the quantum reaction rate to establish the notations used in the rest of this Letter. The starting point is the quantum rate expression derived from correlation functions [ 121: k A_B=k(t)=
30 April1993
real- and imaginary-time propagators in eq. ( 1) are replaced by functional integrals over all paths. Fig. 1 shows a diagrammatic representation of the integrals involved in the calculation of k,,, [ 121.In fig. 1, the hypersurface 8 divides space into A and B states. The solid line represents the action for (41 exp( -BH) 14'), and the dashed and dotted lines the actions for (q’ 1exp(i&/h) (4”) and (4” /exp( i%/A) 1q). The integral runs over all positions q on the A side, all positions q’and q’ on the B side, and finally over all possible paths connecting these three vertices. By a coordinate transformation, one can project out the centroid of the thermal path qo=(@)-’ X16 dtq(t), where q(0) corresponds to q in fig. 1 and q(@i) to q’. VCM factor eq. ( 1) in such a way to focus attention on the role of the thermal path centroid q,,
where p(q)=Q-'l~q(t)d(qo-q)exp{-S[q(t)l} is the equilibrium density for finding the centroid at position q and S is the action for the thermal path, v(q; T) is a dynamical frequency factor which encompassesall the remaining path integralsin eq. ( I ) , and I is defined by the last equality. Apart from a multiplicative constant, the product p(q)v(q; F) is equal to the quantum flux across the dividing surface with the centroid of the thermal path constrained at position q. VCM argue based on analytical considerations that the dynamical frequency factor y( q; f)
-? @ix* Q-’
xImTr[hAexp(-/?H) x h, exp(itH/fi) he exp( -ilH/fi)]
,
(1)
where hA and hg= 1-hA are the population operators of the reactant state and the product state, respectively, H is the Hamiltonian, Q is the canonical partitionfunction,x,=Q-‘Tr[h,exp( -pH)] isthe equilibrium mole fraction of state A and tis the plateau time [ 11. Using Feynman’s path integrals, the
Fig. 1. Diagram representing the path integral in eq. ( 1). The solid line represents an imaginary-time path; the dashed and dotted lines represent real-time paths. The surface @ divides space into A and B states. The integral runs over all positions 4 on the A side, all positions q’and q”on the B side and over all paths.
131
Volume 206, number I ,2,3,4
CHEMICAL.PHYSICSLETTERS
should be sharply peaked at @ such that the dominant contribution to the total flux will come from q near the dividing surface, and the factorization iu the second equality in eq. (2) has been made to emphasize this.
3. Hamiltonian model for electron transfer reactions In the following, we study the nature of the centroid density p(q) and the dynamical frequency factor u(q; f) of electron transfer reactions in the tightbinding limit within a class of simple Hamiltonian models, H=H,
tHB(X) tp[ B(X)t Jo] )
1 dtC(t), --m
W
where C(t)=Q-‘Tr{exp[ xexp( -itH,/h)}
-(fi-it)H,/h] ,
(4b)
and HA and Hu are the diabatic (K=O) Hamiltonians of states A and B. Wolynes [ 91 suggestedthat the integral in eq. (4a) could be more efficiently computed by distorting the integration contour to go from t* -co to r? t co such that the contour will pass through a stationary point t* in the complex-t plane. For the Hamiltonian in eq. (3) with &,=O, a sta132
tionary point does exist at t* = - f ifi, and as has been pointed out in refs. [7,10,11], the value of the integraud C( t* ) is precisely the centroid density p ($ ) on the dividing surface. Therefore, the centroid formulation in the Golden Rule limit appears to correspond exactly to the saddle point approach of Wolynes. Notice, however, that the reaction coordinate considered in these previous studies aud the present one is a discrete variable in the tight-binding limit, whereas the VCM centroid formulation was originally proposed for a continuous reaction coordinate.
4. Quantum dynamics simulations of electron transfer rate
(3)
where Hm= -K&; and & is the Pauli spin operator, HB is a bath Hamiltonian, ,L=,ut$ is the dipole operator, b= Ci cix, is the solvent polarization with xi being the coordinate of the ith bath mode, and &#O for asymmetric electron transfer. In the following, we will treat the symmetriccase &o=0 only. We also assume HB is harmonic, in which case the influence of the bath on the two-level system is captured entirely by the bath spectral density J(o) [ 13,141, which is related to the bath correlation function ( B( 0) J(t) ) [ lo]. Notice that within this tight-binding model, the position of the electron is described by a discrete spin variable 0 - the electron is in the A state if 6= + 1 and in the B state if cr= - 1. In the limit K-+0, the rate can be calculated using the Golden Rule of perturbation theory [ 9,101: xAk~~B=(K/fi)z
30 April 1993
To understand the factorization made in the centroid formulation when it is applied to the tightbinding electron transfer model, we study in great detail the behavior of the centroid density p(q) of the electronic’coordinate and the dynamical frequency factor v(q; t) as a function of q with rigorous quantum Monte Carlo methods. We adopt the aqueous Fe2+/Fe3+ self-exchange reaction as a model, using the Hamiltonian in eq. (3) with the spectral density J(o) derived from the molecular dynamics simulationsof Kuharski et al. [ 15] at room temperature. With the classical isomorphism [ 161, the quantum system is mapped onto a 1D king model which has nonlocal influence functional bonds [ 16181. The calculations of p(q) involve purely imaginary-time path integrals and were carried out using standard umbrella sampling methods. The real-time quantum Monte Carlo simulations of v(q; t ) were performed using a filtering scheme proposed by Mak [ 191 in which the natural weight of each path q(t) in state space is filtered by a damping function D[q(t)]=l+{(1)[q(t)]+<(2)[q(t)], where {(I) [q(t) ] and c(*) [q(t) ] are the sum of the weights of all first-order and second-order fluctuations from the path q(t) defined in ref. [ 191 #I. xl Since the centroid of the imaginary-time path is constrained, not all first- and second-order fluctuations are allowed. The results reported here were obtained using all first- and second-order fluctuations and are therefore subject to a small uncertainty in the centroid position. However, the error introduced is small as long as the number of discretizations in the imaginary-time path is large,andcan bc further reducedby increasingthis number.
CHEMICALPHYSICSLETTERS
Volume 206, number 1,2,3,4
-1
-a.5
0
0.5
1
0
Q
0.2
0.4
30 April 1993
-1
-0.5
Kt
0
0.5
1
9
Fig.2. MonteCarlo results for,W=0.4. (a) Centroid freeenergycomputedusing 100 discretizations for the thermal path. (b) Dynamical frequency factor v(q;t)for several values of q.From top to bottom, q= 0,O.1250.25 and 0.375. 32 discretizations were used for the thermal path and each of the real-time paths. v(q; t) reaches a plateau after a short initial transient. (c) The plateau value v(q; r) as a function of q.
We first examine the results for small electronic coupling jlK= 0.4, which is within the Golden Rule limit. Fig. 2a shows the centroid free energy BF( q) = - In p ( q). As expected, the centroid free energy has a maximum at the dividing surface q= @ = 0, so the centroid density p(q) would have a pronounced minimum at @*.The dynamical frequency factors v(q; t) obtained from quantum dynamics simulations are shown in fig. 2b for several values of q. After an initial transient, v (q; I) for all q approach their plateau values v(q; L). Fig. 2c shows that the plateau values v( 4; t) decay rapidly as a function of q away from the maximum at g=@, and as VCM suggest, v (q; f) is indeed sharply peaked at 8. According to eq. (2), the rate does not depend on either p( q) or v( q; f) individually, but on the product p(q)v(q; f), which is proportional to the quantum flux with the thermal path centroid constrained at q. VCM suggest that because the dynamical frequency factor is sharply peaked at @, the flux should be sharply peaked there as well. As VCM demonstrate, this behavior is certainly true for classical or nearly classical systems [ 71. These are systems with continuous potentials along the reaction coordinate and above the crossover temperature T,: defined by AcoJkT, < 214 where U, is the barrier curvature. But for the highly quanta1 electron transfer system studied herein, the numerical results derived from the data in figs. 2a and 2c and plotted in f%. 3 show a different behavior, They show that the flux is constant within statistical errors for a broad range of q. The sharp decrease of Y(q; f) away from #’ appears to be delicately balanced by the rapidly increasing
p(q) so that the product is invariant with q. Because of the rapidly decreasing magnitude for v( q; f) away from q* = 0, substantially longer sampling times are required to compute v(q; f) with the same relative error. The time required to obtain the last point in fig. 3, which corresponds to q= to.4375 was a factor of 100 longer than q=O, and to further reduce the error to the same level as q=O would require a factor of 2500 more computational effort. To extend our calculations beyond the present limit would require impractically long sampling times. The results in fig. 3 strongly suggest that in the Golden Rule limit, p( q)v( q; f) is invariant with q in the tight-binding limit and is not dominated by the flux at a single point q*, in contrast with VCM’s observation for a continuous reaction coordinate. However, VCM’s suggestion that the rate is ultimately governed by the value of p( q) at q= fl is still
-I
-0.5
0
0.5
I
9
Fig. 3. The quantum fluxp(q)u(q; T) normalized top(q)v(q;i) plotted versus the centroid position q for /X=0.4. Vertical bars indicate one standard deviation. Dashed line denotes Golden Rule limit.
133
Volume206,number1,2,3,4
CHEMlCAL PHYSICS LETTERS
correct in the tight-binding model, at least in the Golden Rule limit. This is because if p( q) v (q; fl is indeed constant throughout the whole range of g, then according to eq. ( 2 ) /c*_~ is simply p (8) multiplied by a constant f I where r= 4 v(8; f). In this case, a short real-time quantum Monte Carlo simulation used to compute v(q; F) at a single point q= 8 plus a straightforwardimaginary-timecalculationfor p( q) are sufficient to compute an exact quantum rate. Having demonstrated numerically that in the Golden Rule limit p(q)v(q; f) is constant, we will now give a mathematical proof for it. The Golden Rule expressionin eq. (4) can be derived easily from eq. ( 1) by retaining only those diagrams that have exactly two kinks - one in the imaginary-time path and the other in one of the real-time paths [ 121. The kink in the imaginary-time path may occur at an imaginary time t anywhere between 0 and m so we can classify contributions to the Golden Rule rate constant according to the location of the kink along the imaginary-time path, (sa) where k”‘(~)=(K/h) xTr{exp[ -
7 dtQ-’ -co (@-7)HAlfi]
xexp( -it&/fi)
exp(itWh)
exp( -.r&)}
.
(5b)
Analytically continuing the integrand into the complex-t plane and using Cauchy’s theorem, it is easily shown that kc’) (7) is invariant with z. Within the two-state model defined in eq. (3)) z is related to the centroid of the imaginary-time path q. by a simple relationship #2:qo= 1- 2~/Bjl.It follows that k(‘) (r) is identical to p( qo)v ( qo; f) apart from a multiplicative constant, and since kc2’(z) is z invariant, p(qo)v ( qo; r) is also invariant with qo. These arguments also establish a clear relationship between a1This is because the time 7 separates the imaginary-time path into two segments, one having all down spins and the other hav-
ing all up spins. The length ofthe first segment is c and the second is J%I- 7. Hence the centroid of the whole imaginary-time path is so=(P)-‘I-r+(Bh-r)l.
134
30April 1993
Wolynes’imaginary-time path integral approach [ 91 and the centroid formulation in refs. [7,10,11]. Having demonstrated that the flux is invariant with the position of the centroid in the Golden Rule limit, we now return to the quantum Monte Carlo results and examine the behavior of p (q) and v(q; f) outside the Golden Rule limit. Fig. 4a shows Monte Carlo results for Y(q; t ) as a function of t for a larger electronic coupling j?K=2, which is well outside the Golden Rule regime. These results and data from centroid density calculations are summarized in table 1. The centroid density p (q) = exp [ - BF( q) ] and the dynamical factor v(q; f) display behaviors similar to the data shown earlier in fig. 2 for /X=0.4. Within statistical error, the quantum flux p( q)v( q) for /X=2 is approximately invariant with the centroidpositionqforO
0
0 0
i Kt
2
0
L
2 Xl
Fig. 4. Monte Carlo results for the dynamical frequency factor: (a) j3K=2, q=O, 0.125,0.25 and 0.375 from top to bottom. 32 discretizations were used for the thermal path and each of the real-time paths; (b) /WA 5, q=O, 0.125 andO.25 from top tobottom. 48 discretizations were used for the thermal path and 24 for each of the real-time paths.
Volume 206,number1,2,3,4
CHEMICALPHYSICS LETTERS
30 April 1993
Table 1 Results of Monte Carlo simulations for flK= 2
ww
4
4%
0 0.125 0.25 0.375
1 OS967 0.1516 0.01323
T)
P(q)
P(f7MK i)/P(W(o; i)
Error (%)
MC passes
0 0.4546 1.8314 4.1516
1 0.9401 0.9464 0.8406
0.71 0.69 2.82 12.5
0.5x 106 1.5 2.5 24.0
exhibits the usual plateau after an initial transient. However, for 4 close to the dividing surface, v (q; t ) displays curious oscillations which appear to persist for rather long times. These oscillations are the signature of coherent quantum processes. Since K determines the intrinsic tunneling frequency of the charge, if K is much larger than the typical bath frequency, the charge cannot be fully solvated, and coherence can persist for an extended duration. But interestingly, the results in fig. 4b also suggestthat the coherence appears to be damped rapidly as the centroid moves away from the dividing surface towards the stable states. The appearance of coherence for large values of K creates serious problems for the direct application of the centroid formulation to simulation of quantum electron transfer rates in the crossover region between the diabatic and adiabatic limits. Since largeamplitude coherent oscillations in the flux may persist for an extended period of time, a long real-time path integral calculation may be necessary to reach the plateau time, and as a result, such quantum simulations may be impractical. At the same time, the appearance of coherent motion implies the prevalence of kink formation within the real-time paths. Therefore, it is not clear whether a method based solely on imaginary-time calculations alone can reproduce an electron transfer rate which is possibly highly influenced by coherent processes. Our results suggestthat the centroid breakup of the full electron transfer rate indicated by eq. (2) in terms of the electronic centroid is perhaps not optimal for systems approaching the adiabatic limit. Based on the physical situation, it may be argued that a similar breakup of the electron transfer rate in terms of the nuclear centroid is more appropriate [20], but the validity of this alternative formulation remains to be tested [2 11. We should point out that although our simulations suggestthat the centroid formulation for
the electron transfer rate using the electronic centroid may not be optimal for systems in between the diabatic and adiabatic limits, a recent calculation [ 221 of the electron transfer rate based on the electronic centroid free energy alone yields an activation energy that appears to be consistent with the conventional estimates in both the diabatic and adiabatic limits [24]. Finally, we also note that for strictly adiabatic systems, a computational scheme based on an entirely classical treatment of the solvent modes is sufficient for computing the electron transfer rate constant [ 231 #3.
Acknowledgement Throughout this work, the authors have benefited from many enlightening discussions with David Chandler, Masako Takasu, Reinhold Egger,Gregory Voth and Joel Bader. CHM wishes to thank David Chandler for introducing him to this area and for suggestingthe studies reported here. We also thank David Chandler and Greg Voth for helpful comments on the manuscript. This research was supported by the Camille and Henry Dreyfus Foundation, the USC FRIF Program and the Army Research Office. Computational resources from IBM are gratefully acknowledged. (iaFor a recent application to an electron transfer reaction in the strictly adiabatic limit, see ref. [ 241. For a related treatment of adiabatic proton transfer, see ref. [ 251.
References [ 1I D.Chandler, J. Chem. Phys. 68 ( 1978) 2959. [2 J J.T. Hynes,irx The thcoq of chemicalreactions,ed. M. Baer (CRC Press, Bcca Raton, 1985).
135
Volume 206, number 1,2,3,4
CHEMICALPHYSICSLETTERS
[ 31 B.J. Beme, in: Multiple time scales, eds. J.V. Brackbilland 9.1. Cohen (Academic Press, New York, 1985) . [ 41 R.O. Rosenberg, B.J. Beme and D. Chandler, Chem. Phys. Letters 75 (1980) 162. [ 5] C.H. Bennett,in: Algorithmsfor chemicalcomputation, ACS Symp.Ser. No. 46, ed. R.E. Cbristofferson (Am. Chem. Sot., Washington, 1977). [6] T. Yamamoto, J. Chem. Phys. 33 (1960) 281. [ 71 G.A. Voth, D. Chandler and W.H. Miller, J. Chem. Phys. 91 (1989) 7749. [8] M.J. Gillan, J. Phys. C 20 (1987) 3621. [ 9) P.G. Wolynes,J. Chem. Phys.,87 ( 1987) 6559. [lo] D. Chandler, in: Liquids, freezing and the glass transition, Les Houches 51, Part 1, eds.D. Levesque, J.P. Hansen and J. Zinn-Justin (Elsevier, Amsterdam, 1991). [ 111J.S. Bader, R.A. Kuharshi and D. Chandler, J. Chem. Phys. 93 ( 1990) 230. [ 121G.A. Voth, D. Chandler and W.H. Miller, J. Phys. Chem. 93 (1989) 7009. [ 131A.O.Caldeiraand A.J. Leggett,Phys. Rev. Letters 46 ( 1931) 211.
136
30 April 1993
[ 141A.J. Leggett, S. Chakravarty, A.T. Dorsey, M.P.A. Fischer, A. Garg and W. Zwerger, Rev. Mod. Phys. 59 (1987) 1. [ 151R.A. Kuharski,J.S. Bader,C. Chandler,M. Sprik,ML. Klein and R.W. Impey, J. Chem. Phys. 89 (1988) 3248. [ 161D. Chandler and P.G. Wolynes,J. Chem. Phys. 74 (1981) 4078. [ 171M. Suzuki,Pros. Theoret. Phys. 56 (1976) 1454. [ 181D. Chandler, Introduction to modem statistical mechanics (Oxford Univ. Press, Oxford, 1987) [ 191C.H. Mak, Phys. Rev. Letters 68 (1992) 899. [ 201D. Chandler, private communication. [ 211C.H. Mak, in preparation. [ 22 1J.N. Gehlen, D. Chandler,H.J. Kim and J.T. Hynes, J. Phys. Chem. 96 ( 1992) 1748. [ 231R.A. Marcus and N. Sutin, Biochim. Biophys. Acta 8 11 (1985)265. [ 241D.A. Zichi, G. Ciccotti, J.T. Hynes and M. Ferrario, J. Phys. Chem. 93 (1989) 6261. [ZS] D. Borgis, G. Tarjus and H. Azzouz, J. Phys. Chem. 96 (1992) 3188.