J. Phys. Chem. 1996, 100, 1137-1143
1137
Properties of Quantum Transition State Theory and Its Corrections Nancy Fisher Hansen and Hans C. Andersen* Department of Chemistry, Stanford UniVersity, Stanford, California 94305 ReceiVed: June 12, 1995; In Final Form: September 22, 1995X
A previously proposed version of quantum transition state theory (QTST) is analyzed in the classical (p f 0) limit. In this limit, the theory becomes classical transition state theory for the particular transition state (TS) used. In addition, a variational choice for the transition state surface for use in this theory is suggested. An alternative quantum transition state theory (QTST2) is presented, and results are presented for symmetric and asymmetric Eckart barriers. Corrections to both theories based on analytic continuation are proposed. QAk ) ∫0 Cf(t) dt ∞
I. Introduction Recently we proposed a new quantum mechanical transition state theory (QTST) for adiabatic reactions.1 On the basis of the flux-flux correlation function of Miller, Schwartz, and Tromp,7 QTST predicts rates by integration of the exact correlation function for a one-dimensional parabolic barrier with frequency and height chosen to match the zero time value and second derivative of the correlation function calculated for the true potential. It is one of many proposed versions of quantum transition state theory.2-6 The purpose of this paper is to examine various properties of this theory and expand it to include a prescription for choosing a transition state (TS) surface. We will show that in the classical (p f 0) limit, our theory becomes classical transition state theory8 for the particular transition state surface used. The proposed choice of TS surface is based on an analysis of the parabolic barrier linearly coupled to a harmonic oscillator. In the classical limit, this surface becomes the one that minimizes the predicted value of the rate; that is, it becomes the optimal transition state surface for use in classical transition state theory. Another quantum transition state theory (QTST2), similar to QTST in its structure and simplicity, is also presented. On symmetric and asymmetric Eckart barriers,9 QTST2 gives answers similar to QTST. In particular, it overestimates the rate constant for highly asymmetric barriers at low temperatures. Various correction schemes based on analytic continuation schemes are applied to both theories. Section II contains a derivation of the classical limit of QTST. Section III examines the dependence of QTST on the choice of transition state surface, and an optimal surface is proposed. In section IV, the new, second version of quantum transition state theory is described, and results are presented for asymmetric Eckart barriers. Section V is a discussion of corrections to both theories, and a concluding discussion is in section VI. II. Classical Limit of QTST A. Statement of QTST. For a system with reaction coordinate Q and Hamiltonian
H)
P2 p2 + + V(Q,q) 2m 2m
(1)
(where q represents any additional degrees of freedom and p represents the conjugate momenta), the exact rate constant k for the conversion of reactants to products is given by X
Abstract published in AdVance ACS Abstracts, January 1, 1996.
0022-3654/96/20100-1137$12.00/0
(2)
where QAk is the partition function for reactant states and Cf(t) is the flux-flux correlation function,7
Cf(t) ) Tr[e-βH/2Fe-βH/2F(t)]
(3)
Here β ) 1/kBT, F is the flux operator,
F)
P P 1 δ(Q - Q†) + δ(Q - Q†) 2 m m
(
)
(4)
Q is the reaction coordinate, and P is its conjugate momentum. Here it is assumed that the surface defined by Q ) Q† is the dividing surface in configuration space between reactant and product states, reactants having Q < Q† and products having Q > Q†. For a particle moving in one dimension across a parabolic barrier Vpar(s) ) Eb - 1/2mωb2s2, the flux-flux correlation function is7
Cf,par(t;ωb,Eb) )
ωb2 sin(βpωb/2) cosh(ωbt)e-βEb 4π(sin2(βpωb/2) + sinh2(ωbt))3/2
(5)
To obtain a theory that is exact for a one-dimensional parabolic barrier, it is advantageous to construct the function
R(t;ωb,Eb) ≡
Cf(t) Cf,par(t;ωb,Eb)
(6)
Here, ωb and Eb are the frequency and barrier height for a onedimensional parabolic reference system chosen for the particular problem at hand. Then the rate constant is given by
QAk ) ∫0 Cf,par(t;ωb,Eb) R(t;ωb,Eb) dt ∞
(7)
Equation 7 is exact for any choice of ωb and Eb, but to facilitate the analytic continuation of R, ωb and Eb can be chosen to make R as simple and as slowly varying as possible. It is convenient to choose them such that
R(t;ωb,Eb) ) 1 + O(t4) This is equivalent to demanding that © 1996 American Chemical Society
(8)
1138 J. Phys. Chem., Vol. 100, No. 4, 1996
Hansen and Andersen
Cf,par(0;ωb,Eb) ) Cf(0)
(9)
C′′f,par(0;ωb,Eb) ) C′′f(0)
(10)
That is, the reference system must have the same value and the same second derivative of Cf at t ) 0 as the system of interest. (The odd derivatives of both functions at t ) 0 are zero.) The correlation function and its second derivative at t ) 0 can be expressed in terms of values of the density matrix and its derivatives with respect to the coordinates, evaluated on the transition state surface. QTST consists of making the simplest approximation for R(t) by assuming that R(t) ) 1 for all t. The corresponding expression for the rate is
C ¨ f(0) )
()
(
where V(n,m)(Q†,q) denotes the mixed partial derivative ∂n+mV(Q,q)/ ∂Qn∂qm evaluated at the point (Q†,q). In the classical limit, the exact rate constant of a parabolic barrier becomes
ZAkQTST )
QAkQTST ≡ ∫0 Cf,par(t;ωb,Eb) dt ) (QAk)par,ωb,Eb (11)
Then the exact rate can be expressed as
QAk ) QAkQTSTκ
(12)
)
where
κ ) ∫0 Cf,par(t;ωb,Eb) R(t;ωb,Eb) dt/∫0 Cf,par(t;ωb,Eb) dt
() 1
β2p2ω2b
is the dynamical correction factor or transmission coefficient for this theory. B. Classical Limit Derivation. We now show that this formulation of QTST, for any choice of dividing surface, gives a rate constant that in the classical limit (p f 0) approaches the classical transition state theory8 rate constant for the same choice of dividing surface. An analysis of this theory in the classical limit involves ¨ f(0) in this limit. Makri determining the behavior of Cf(t) and C and Miller10 presented a series representation of Cf(t) in the limit of small t and small p for a system with F degrees of freedom:
exp{(i/p)W(Q,q,Q′,q′;tc)} 〈Qq| exp(-iHtc/p)|Q′q′〉 ) t-F/2 c (14) where
W(Q,q,Q′,q′;tc) ) ∑tn-1 c Wn(Q,q,Q′,q′)
(15)
n)0
This representation can be used to derive a series expansion of ¨ f(0) in powers of p, which yields Cf(0) and C
{ (
)}
p2 1 1 F-1 dq dp exp -β + V(Q†,q) × ∫ 2 2 h 2m πβ p β2p2 (2,0) † β2p2 (0,2) † 1V (Q ,q) V (Q ,q) + 12m 12m β3p2 (1,0) † V (Q ,q) V(1,0)(Q†,q) + 96m β3p2 (0,1) † V (Q ,q) V(0,1)(Q†,q) + O(p4) (16) 24m
(
{(
F-1
∫dq dp exp
)
pi2
)}
-β V(Q ,q) + ∑ i 2m †
×
(1 + O(p2)) (19)
∞
(13)
and
Cf(0)4π sin2(βpωb/2)
h
()
(18)
e-βEb ≡ Cf(0)/Cf,par(0;ωb,0) )
Cf(0) )
kBT -βEb (1 + O(p2)) e h
A consequence of eqs 9 and 10 is that
kBT βpωb e-βEb h 2 sin(βpωb/2)
∞
)} )
∞
)
{ (
p2 -12 1 F-1 dq dp exp -β + V(Q†,q) × ∫ 2m πβ4p4 h β2p2 (2,0) † β2p2 (0,2) † 1V (Q ,q) V (Q ,q) + 12m 9m 7β3p2 (1,0) † V (Q ,q) V(1,0)(Q†,q) + 96m β3p2 (0,1) † V (Q ,q) V(0,1)(Q†,q) + O(p4) (17) 12m
Therefore, the classical limit of our quantum transition state theory is
lim ZAkQTST ) pf0
()
kBT 1 h h
F-1
{(
∫dq dp exp
pi2
)}
-β V(Q ,q) + ∑ (1 + O(p2)) i 2m (20) †
This is simply the classical transition state theory rate constant for the dividing surface Q ) Q†. The value of ωb is irrelevant in the classical limit. III. Variation of the Transition State Surface Classical transition state theory is guaranteed to predict a value higher than the true classical rate constant. For this reason, the best choice for a transition state surface for use with this theory is obviously the one which predicts the lowest rate constant. Quantum transition state theory, however, does not produce an upper bound to the true quantum rate constant,1 and therefore the best choice of a transition state is not as clear. In many cases, it can be indicated by symmetry considerations, but in many others, there are a number of surfaces that could reasonably be chosen. We examine the parabolic barrier linearly coupled to a harmonic oscillator to see the effect of changing the TS surface, which in this case is a line in (Q,q) space. The potential for this system is
1 1 V(Q,q) ) - MΩ2Q2 + mω2q2 + CqQ 2 2
(21)
We now want to consider the predicted rate constant for this potential when the dividing surface is an arbitrary straight line, Q ) Q† + φq. Each choice of a potential of this type (specified by the dimensionless parameters R ) βpΩ, γ ) C/(Mm)1/2ωΩ, and δ ) ω/Ω) and a TS of the form Q ) Q† + φq can be
Properties of Quantum Transition State Theory
J. Phys. Chem., Vol. 100, No. 4, 1996 1139 singularities of the appropriate type at the appropriate places on the imaginary axis. QTST2 is based on such a construction. To simplify matters, let
C ˆ f(z) ) Cf(z1/2)
(24)
ˆ f(t2) Cf(t) ) C
(25)
or Figure 1. Region of definition of Cf(x), within which Cf is necessarily analytic. For a parabolic barrier with ωb > π/βp, however, there are singularities within this region, at t ) (π/2ωb.
changed by a point transformation into an uncoupled parabolic barrier and a TS of the same form. Varying the TS in the new uncoupled system is equivalent to varying it in the old and is much simpler. So we now consider potentials of the form
V(Q,q) ) DQ2 + Eq2
(22)
where D < 0 and E > 0, with dividing surface
Q ) Q + Φq †
Any approximate representation of C ˆ f(z) will automatically lead to a representation of Cf(t) that is even in t, as is true for any physical Cf. The discontinuous population operators in the flux operator imply the existence of a branch point of C ˆ f(z) at z ) ˆ f in the vicinity of these -β2p2/4 and the following behavior of C points.11
C ˆ f(z) ≈
βpλ〈Q†| exp(-βH)|Q†〉
+ O(z + β2p2/4)-1/2 (26)
8π(z + β p /4) 2 2
3/2
Now consider the function
(23)
¨ f(0), and ΓQTST for these systems. and calculate Cf(0), C For a parabolic barrier with one uncoupled harmonic oscillator, QTST with dividing surface Q ) 0 gives the exact rate. This choice of TS corresponds to a saddle point in the predicted rate constant with respect to variation of Q† and Φ. The QTST rate constant is a minimum with respect to moving the TS parallel to the reaction coordinate, and a maximum with respect to rotation of the TS. Cf(0), however, is a minimum with respect to rotation or translation of a linear TS from the one which gives the exact rate constant. In addition, Cf(0) determines the value of the rate constant in the classical (p f 0) limit, so the best choice of a linear TS for classical transition state theory is the location of a minimum of Cf(z). Motivated by these facts, we suggest that the best choice of a dividing surface for QTST is the one that minimizes Cf(0). When this choice is made, QTST is exact not only for the onedimensional parabolic barrier, but also for any system of a parabolic barrier linearly coupled to harmonic oscillators. IV. An Alternative Version of Quantum Transition State Theory A. Basic Theory. QTST makes use of the value and second derivative of Cf(t) at t ) 0 to estimate the integral of the function from zero to infinity and hence to estimate the rate constant. The calculation makes use of the exactly known behavior of Cf,par for the parabolic barrier problem. Although it is not necessary to describe the theory in these terms, QTST can be regarded as, in effect, making the approximation that Cf(t) ) Cf,par(t) for all t, for some suitably chosen parabolic potential. When viewed in this way, QTST has a troubling aspect, namely, that the resulting approximation for Cf(t), regarded as a function of complex t, has branch points within the region where it should be analytic (see Figure 1). This occurs when the frequency ωb of the parabolic reference system is greater than π/βp. The branch point is an artifact of the fact that the parabolic barrier is un unphysical mechanical system, with a potential energy function that is not bounded from below and that approaches -∞ as the coordinate goes to (∞. As an alternative to QTST, we might imagine using the same information about the true correlation function at t ) 0 to construct an alternative approximation formula for Cf(t) that not only has all the symmetry properties of the true Cf but that also is analytic in the appropriate region of analyticity and has
D(z) ≡
d ln(C ˆ f(z)) dz
(27)
At the physical singularities of C ˆ f, D has poles with residue -3/2. In addition, D will have poles wherever C ˆ f has zeros. If we assume that C ˆ f has only the physical singularity and no zeros (this is the case for the free particle as well as the parabolic barrier problem), then a reasonable Pade´ approximant for D would be
D(z) ≈
a′ + b′z 1 + 4z/(β2p2)
(28)
where
b′ )
24 + 6a′ β2p2
(29)
so the residue will be correct. Integrating D(z) to obtain a model C ˆ f function yields
C ˆ f,mod(z;a,b) )
a exp(bz)
(30)
(z + β2p2/4)3/2
(The constant of integration has been combined with a′ to give the new constant a and b ) 6 + a′). The approximation made in this alternate version has many desirable characteristics. Provided b is negative, it will always be integrable and give a positive rate constant. It has a branch point in the right location and is analytic where we know C ˆ f should be. Now, in the spirit of QTST, we construct the function
R2(z;a,b) )
C ˆ f(z)
(31)
C ˆ f,mod(z;a,b)
As before, the exact rate constant is given by
ˆ f,mod(z;a,b) ∞C
QAk ) ∫0
R2(z;a,b)
1/2
2z
dz
(32)
In the same way we did with QTST, we choose the parameters a and b so that near z ) 0,
R2(z;a,b) ) 1 + O(z2)
(33)
that is, the C ˆ f,mod should have the same value and derivative at
1140 J. Phys. Chem., Vol. 100, No. 4, 1996
Hansen and Andersen
TABLE 1: QTST2 Tunneling Correction Factor Γ for Symmetric Eckart Barriersa ) 24
) 48
) 72
ξ
Γexact
ΓQTST
ΓQTST2
Γexact
ΓQTST
ΓQTST2
Γexact
ΓQTST
ΓQTST2
2 4 6 8 10 12
1.22 2.04 3.69 6.15 9.05 12.1
1.24 2.12 3.89 ... ... ...
1.29 2.20 4.03 ... ... ...
1.21 2.07 4.54 10.5 22.5 42.4
1.22 2.21 5.05 11.4 23.3 41.9
1.26 2.21 4.92 11.2 23.0 42
1.20 2.10 5.20 15.6 45.7 118
1.20 2.25 6.08 18.4 51.3 124
1.25 2.19 5.61 16.8 47.2 116
ξ
Γexact
ΓQTST
ΓQTST2
Γexact
ΓQTST
ΓQTST2
Γexact
ΓQTST
ΓQTST2
2 4 6 8 10 12
1.20 2.11 5.75 21.8 87.3 307
1.20 2.26 6.93 27.2 104 340
1.25 2.16 6.16 23.7 90.7 301
1.20 2.13 6.23 29.4 162 783
1.20 2.27 7.66 38.1 199 888
1.24 2.14 6.63 32.1 168 762
1.19 2.14 6.66 38.8 295 1970
1.21 2.28 8.29 51.6 369 2269
1.24 2.12 7.03 42.5 304 1900
) 96
a
) 120
) 144
Exact results are from Johnston, ref 9. (...) indicates that a positive value of b prevented integration of C ˆ f.
the origin as does the true C ˆ f. This can be accomplished by requiring
a)
b)
βpC ˆ f(0) 8
)
βpCf(0) 8
ˆ f(0)/β2p2 C ˆ ′f(0) + 6C C ˆ f(0)
)
C ¨ f(0) 6 + 2Cf(0) βp 2 2
(34)
The assumption of this alternate version of quantum transition state theory (QTST2) is that R2(z;a,b) ) 1 for all z. The QTST2 expression for the rate is then
ˆ f,mod(z;a,b) ∞C
QAkQTST2 ) ∫0
2z1/2
dz
(35)
where a and b are defined by eq 34. The dynamical correction factor for this theory is
ˆ f,mod(z;a,b) ∞C
κQTST2 ) ∫0
R2(z;a,b)
2z1/2
ˆ f,mod(z;a,b) ∞C
dz/∫0
2z1/2
B. Classical Limit of QTST2. An analysis of QTST2 in the classical (p f 0) limit reveals that, like QTST, it becomes classical transition state theory for the dividing surface used. In this limit,
b ) O(p2)
(37)
a ) βpCf(0)/8
(38)
ˆ f: and C ˆ f,mod becomes the free particle C
Cf(0) 8(z + 1/4)3/2
(1 + O(p2))
barriers with QTST and exact results is contained in Table 1. For these potentials,
dz (36)
C ˆ f,mod(z) )
Figure 2. Modified Arrhenius plot of the QTST2 rate constant for the parabolic barrier. The vertical axis is ln(βhZAk), and the horizontal axis is the barrier height divided by kBT. The solid curve is the exact results, the dashed curve is the QTST2 prediction, and the dotted curve is the prediction of classical transition state theory.
(39)
The integral of this, as in eq 39, gives the rate constant of classical transition state theory. Unlike QTST, QTST2 does not predict the exact correction factor for a one-dimensional parabolic barrier. It does, however, predict the divergence at ξ ) βpω ) 2π. As can be seen in Figure 2, QTST2 overestimates the rate constant at higher temperatures and underestimates them when the value of ξ is larger. C. Results for Symmetric and Asymmetric Eckart Barriers. A comparison of the QTST2 results for symmetric
V(Q) ) V0 sech2(Q/∆)
(40)
The tunneling correction Γ,12 defined by
Γ)
(ZAk)quantum (ZAk)classical
) 2πβp(ZAk)quantum
(41)
is a function of two dimensionless parameters, ) βV0 and ξ ) βp(2V0/m∆2)1/2, where m is the mass of the particle. The density matrix was calculated using the numerical matrix multiplication algorithm of Thirumalai, Bruskin, and Berne.13 The Arrhenius plot is given in Figure 3. For many of the symmetric Eckart barriers considered, QTST2 predicts a rate constant slightly better than QTST, and in all cases it makes a comparable prediction. We now consider an asymmetric Eckart barrier
V(Q) )
(V1 - V2) 1 + exp(-2πQ/∆)
+
(V11/2 + V21/2)2 4 cosh2(πQ/∆)
(42)
An asymmetric Eckart barrier is completely specified by V1, the height from the reactant region to the top of the barrier, V2, the height from product region to the barrier, and ∆, which is approximately the width. We will consider different temperatures with the dimensionless parameter
Properties of Quantum Transition State Theory
J. Phys. Chem., Vol. 100, No. 4, 1996 1141
Figure 3. Modified Arrhenius plot of the rate constant for the symmetric Eckart barrier with R(≡2π/ξ) ) 12. The vertical axis is ln(βhZAk), and the horizontal axis is the barrier height divided by kBT. See the text for the meaning of the symbols. The squares are the QTST results, and the diamonds, which closely coincide with the squares, are the QTST2 results. The solid curve is the exact result. The dashed straight line gives the classical transition state theory result. The QTST2 rate constants, like the QTST and exact results, display a gradual transition from activated behavior at high temperatures to a less steep temperature dependence, characteristic of tunneling, at lower temperatures.
Figure 5. Modified Arrhenius plot of the rate constant for the asymmetric Eckart barrier with R1 ) 12, R2 ) 48. The meaning of symbols is the same as in Figure 4.
ZAk ) ∫0 dt Cf,app(t) R(t)
(44)
R(t) ≡ Cf(t)/Cf,app(t)
(45)
∞
where
2. Choose the parameters of Cf,app so that R will be one with zero slope and second derivative at t ) 0. 3. Make the approximation that R(t) ) 1 for all t, so
ZAkQTST ) ∫0 dt Cf,app(t) ∞
Figure 4. Modified Arrhenius plot of the rate constant for the asymmetric Eckart barrier with R1 ) 12, R2 ) 24. The vertical axis is ln(βhZAk), and the horizontal axis is the barrier height divided by kBT. See the text for the meaning of the symbols. The squares are the QTST results, and the diamonds are the QTST2 results. The solid curve is the exact result. The dashed straight line gives the classical transition state theory results.
ξ)
(
)
2 βp 8π V1V2 m1/2 ∆2
1/2
1 (V11/2 + V21/2)
(43)
As for the symmetric Eckart barriers, we will consider constant potential trends specified by the parameters R1 ) 2πβV1/ξ and R2 ) 2πβV2/ξ. The maximum of the potential is no longer located at zero but at Qmax ) ∆ ln(V1/V2)/(4π). As suggested in section III, the transition state Q† is chosen to be the one that minimizes the value of Cf(0). Both QTST and QTST2 do not do well for highly asymmetric Eckart barriers in the low-temperature limit. Figures 4 and 5 show modified Arrhenius plots for Eckart barriers with R1 ) 12, R2 ) 24 and R1 ) 12, R2 ) 48, respectively. The QTST2 results are slightly better than the QTST results in this limit. Still, they are not good, and these asymmetric barriers are where we concentrate our correction efforts for both theories. It should be noticed that there is not a large difference between the values when the TS is chosen to minimize the rate constant or Cf(z). Efforts at corrections by analytic continuation alleviate at least some of the problem (see section V). V. Corrections to QTST Each of the two quantum transition theories (QTST and QTST2) can be formulated as follows: 1. Express the rate constant in the form
(46)
In QTST, the approximate Cf is the correlation function for a one-dimensional parabolic barrier Cf,par, and in QTST2, it is Cf,mod, defined in eq 30 as a function of z ) t2. The obvious step to take in making corrections to either theory is to improve our approximation for R for positive time. We will attempt to do this by analytic continuation. When a correction method applies equally well to both theories, we will use the general notation above. A. Analyticity Properties of R. The evenness of R as a function of t suggests14 that we deal with two new functions ˆ f(z) defined as in section Rˆ (z) defined by Rˆ (z) ≡ R(z1/2) and C IV. In performing analytic continuation of Rˆ , it is worthwhile noting its analyticity properties. Rˆ is analytic wherever Cf and Cf,app are both analytic and Cf,app is not zero. At points where Cf,app ) 0, Rˆ will have a pole. But Rˆ (z) will also be analytic at ˆf the “physical singularity” of C ˆ f, i.e., at z ) -β2p2/4, since C and both model functions behave as
(
C ˆ f,app ) z +
)
β2p2 4
-3/2
A(z)
(47)
where A(z) is an analytic function at z ) -β2p2/4. On the basis of these reasons, we attempt to model Rˆ (z) with a Pade´ approximant. If we can calculate one extra piece of information about the function, appropriate Pade´s would be
Rˆ pade(z) ) Rˆ a(z) ) 1 + az2
(48)
1 1 - az2
(49)
and
Rˆ pade(z) ) Rˆ b(z) )
Each of these functions has potential problems. If a < 0, Rˆ a could potentially predict a negative rate constant. If a > 0, Rˆ b would have poles along the positive axis, which we know is not the case.
1142 J. Phys. Chem., Vol. 100, No. 4, 1996
Hansen and Andersen
We therefore make the following prescription: If our estimate of a is positiVe, we use Rˆ a as the Pade´ for Rˆ , and if a is negatiVe, we use Rˆ b. Another alternative is to model Rˆ (z) by an exponential:
Rˆ exp(z) ) exp(az ) 2
(50)
This method can only be used if the second derivative of Rˆ at z ) 0 is negative and therefore will either lower the predicted QTST rate or give no answer at all. Note that this method is equivalent to using a Pade´ approximant, with a single piece of information, for the analytic continuation of ln Rˆ . A consequence of these choices is that our model Rˆ (z), and the resulting model C ˆ f(z), is positive for all positive, real values of z. This is not always the case for the actual C ˆ f(z), so these assumptions may place a limit on the accuracy of the correction. They do avoid the possibility of predicting negative rates, though. To find a value of a for use in the above correction schemes, we need to calculate some additional information about Rˆ (or equivalently, C ˆ f). We choose to take advantage of our ability to calculate the value of C ˆ f along the negative z axis. (The use of C ˆ f(z) along the negative real axis to predict its behavior along the positive real axis is equivalent to using Cf(t) for imaginary time to predict its behavior for real time by analytic continuation in the complex t plane.) If xz ) iτ and τ is real,
C ˆ f(z) ) Cf(iτ) )
(
∂ p2 2 dq dq′ 〈Qq|e-(β/2+τ/p)H|Q′q′〉 × 2 ∫ ∂Q 4M ∂ 〈Q′q′|e-(β/2-τ/p)H|Qq〉 ∂Q′ ∂2 〈Q′q′|e-(β/2+τ/p)H|Qq〉 ∫dq dq′〈Qq|e-(β/2-τ/p)H|Q′q′〉∂Q∂Q′ 2 ∫dq dq′∂Q∂∂Q′〈Qq|e-(β/2-τ/p)H|Q′q′〉〈Q′q′|e-(β/2+τ/p)H|Qq〉 -
)
(51) The matrix elements and derivatives can be calculated in the same way, and with the same ease, as they are for C ˆ f(0). B. The Second Derivative of R ˆ . We know that Rˆ is of the form
Rˆ (z) ) 1 + O(z2)
(52)
so the most obvious piece of additional information to calculate about Rˆ would be its second derivative at z ) 0. Then, according to the schemes above,
1 d2 Rˆ (z)|z)0 (53) 2dz2 There are a variety of ways to fit the second derivative given values of Rˆ along the negative z axis. The first method we implemented (method I) was to fit 1 + az2 to the value of Rˆ at a single point. The point we chose was z ) -z1/2, where z1/2 is the positive value of z where C ˆ f,app has dropped to half of its value at z ) 0. The rationale for this choice is that the behavior of Rˆ for values of z between about 0 and 2z1/2 is of crucial importance for the value of the integral that gives the rate, and we want the information that we use for the negative z behavior to contain information about Rˆ in a region similarly close to the origin. The second method (method II) was to do a least squares fit of 1 + az2 to data for Rˆ (z) for a range of z values along the a)
TABLE 2: Corrected Quantum Correction Factors for the r ) 12 Symmetric Eckart Barriera ξ
QTST version
2 QTST
QTST result 1.2
QTST2 1.3 6 QTST
6.1
QTST2 5.6 12 QTST
2.3 × 103
QTST2 1.9 × 103
I Γcor,pade I Γcor,exp
II Γcor,pade
III Γcor,pade
II Γcor,exp
III Γcor,exp
1.2 1.2 1.3 1.3 5.3 5.1 5.4 5.3 1.8 × 103 1.7 × 103 1.9 × 103 1.9 × 103
1.2 1.2 1.3 1.3 5.2 5.0 5.4 5.3 1.8 × 103 1.7 × 103 1.9 × 103 1.9 × 103
1.2 1.2 1.3 1.3 5.3 5.2 5.3 5.3 1.8 × 103 1.8 × 103 1.9 × 103 1.9 × 103
Γexact 1.2 1.2 5.2 5.2 2.0 × 103 2.0 × 103
aΓ ´ s in eq 49 for Rˆ (z), cor,pade’s were calculated using one of the Pade and Γcor,exp’s were calculated using the exponential function in eq 50. The superscript (I, II, or III) refers to the method used to calculate the second derivative of Rˆ at z ) 0.
negative real axis. For QTST2 that range was from z ) -β2p2/4 to z ) 0. For QTST the range was for z above the first branch point of C ˆ f,par, which is at zb > -β2p2/4 when βpωref > π. For these cases, only points where z > zb were used. Finally (method III), we did a least squares fit of Rˆ (z) to 1 + az2 + bz3, to see how this changed our value of a. All three methods gave comparable results. The resulting corrected rates are shown in Tables 2 and 3. In almost all cases, both Rˆ pade and Rˆ exp result in a better approximation to the quantum correction factor. Note that all three correction methods give consistently similar results. This illustrates that analytic continuation methods for calculation of the rate constant are very stable at this level of complexity. (However, see below.) Moreover, the methods make significant corrections to the rate only when the corresponding quantum transition state theory is significantly in error. Thus, any of these correction methods appears to provide a diagnostic for determining whether the transition state theory is accurate. Method I is the easiest of the methods to implement, since it requires the evaluation of C ˆ f for only one point on the negative z axis. However, it gives results for the rate constant that are consistently close to those obtained by the more elaborate methods II and III. Moreover, in rate calculations for more complicated systems with many degrees of freedom, it would be very worthwhile if an evaluation of C ˆ f(z) at a single point on the negative z axis could provide not only a corrected rate but also a diagnostic for inaccuracies in the transition state theory. Further corrections of a similar type do not improve the answer further. We have tried various extended correction schemes. In general, our experience is that if the scheme forces the approximate C ˆ f to be nonnegative and nonsingular on the positive z axis, then the calculated rate is not significantly better than the values obtained using the methods of this chapter. On the other hand, if the approximate C ˆ f is allowed to be negative, the calculated rate will often be negative, and if the approximate C ˆ f or Rˆ is fit with a Pade´ approximant that can be singular on the positive real axis, then the fitted function is likely to have such a singularity. These observations highlight the ambiguities and extreme sensitivity associated with attempts to perform analytic continuation by numerical methods. They also suggest that the actual C ˆ f for some problems under highly quantum conditions becomes negative on the positive real axis and that this is the reason why analytic continuation methods that force C ˆ f to remain positive will consistently overestimate the rate.
Properties of Quantum Transition State Theory
J. Phys. Chem., Vol. 100, No. 4, 1996 1143
TABLE 3: Corrected Quantum Correction Factors for the r1 ) 12, r2 ) 48 Asymmetric Eckart Barriera ξ
QTST version
2
QTST
1.2
QTST
1.2
QTST
2.1
QTST2
2.0
QTST
6.3
QTST2
5.2
4
6
8
10
12
QTST result
QTST
51
QTST2
38
QTST
1.5 × 103
QTST2
1.1 × 103
QTST
11.8 × 104
QTST2
8.7 × 104
I Γcor,pade
II Γcor,pade
III Γcor,pade
I Γcor,exp
II Γcor,exp
III Γcor,exp
1.2 1.2 1.2 1.2 2.1 2.1 2.2 ... 5.1 4.9 5.9 ... 36 34 38 38 1.0 × 103 0.9 × 103 1.0 × 103 1.0 × 103 8.3 × 104 7.8 × 104 8.1 × 104 8.0 × 104
1.2 1.2 1.2 1.1 1.9 1.9 2.1 ... 5.1 4.9 5.3 ... 37 34 38 38 1.0 × 103 0.9 × 103 1.1 × 103 1.1 × 103 8.3 × 104 7.8 × 104 8.5 × 104( ) * 8.5 × 104( )
1.2 1.2 1.2 1.2 2.8 ... 2.3 ... 5.8 5.7 5.9 ... 37 34 38 38 1.0 × 103 0.9 × 103 1.0 × 103 1.0 × 103 8.5 × 104 8.0 × 104 8.2 × 104( ) * 8.1 × 104( )
*
Γexact 1.2 1.2 2.0 2.0 5.3 5.3 26 26 0.25 × 103 0.25 × 103 0.41 × 104 0.41 × 104
*
a
See Table 2 for an explanation of the notation. (...) refers to points for which the exponential (eq 50) was unusable due to a positive value of Rˆ ′′(0). Asterisks ( ) indicate that values for R were used only in the region in which they could be accurately obtained, i.e., z > -0.0225.
*
VI. Conclusions Starting with an exact expression for the quantum rate constant of an adiabatic process, we have formulated two versions of quantum transition state theory (QTST and QTST2) that have the following characteristics: (1) The QTST rates are obtained from the exact quantum rate theory in much the same way as the classical transition state theory rate can be obtained from the exact classical theory, namely, by extrapolating the behavior of a certain function at short times to all times. (2) The QTST rates can be obtained from calculations of two static equilibrium averages for the problem, these averages being the value and second derivative of the flux-flux correlation function at zero time. (3) In the classical limit, the QTST rates approach the classical transition state theory rate for the same choice of dividing surface. (4) In some cases, a variational criterion can be used to obtain the dividing surface, but the calculated QTST rate is not an upper bound to the correct result. (5) The QTST gives very accurate rates for certain model problems including symmetric and asymmetric Eckart barriers and parabolic barriers coupled to a single harmonic oscillator. For the latter problem, the QTST rate in particular is generally above the exact result and becomes more accurate as the problem becomes more quantum mechanical. (6) The QTST rate is exact for a onedimensional parabolic barrier problem, by construction, if the dividing surface is chosen to be at the top of the barrier. Moreover, we have presented a method of choosing a dividing surface for use in these theories which makes the first theory exact for a parabolic barrier linearly coupled to a harmonic bath. Finally, we suggest several methods of analytic continuation ˆ f,app(z) to give an approximate quantum of the ratio Rˆ (z) ) C ˆ f(z)/C dynamical correction factor for the theory. These correction methods avoid the problems associated with calculating quantum correlation functions for real time (e.g., the “sign problem” associated with path integral dynamics calculations) but for the problems we tested almost always give an improvement in the rate, especially for problems where the transition state theory does badly. Furthermore, all of them give similar predictions for the corrected rate, leading us to believe that they are a stable
indicator of where our theory is in error. Since in our methods we have made the assumption that C ˆ f(z) is positive along the positive real z axis, we suspect that the logical next step would be to amend our method to allow C ˆ f to become negative along the positive z axis while still guaranteeing a positive rate constant. Acknowledgment. This research was supported by National Science Foundation Grant CHE-89-18841 and CHE-94-21884 and benefited from computational resources provided by NSF Grant CHE-88-21737. The graduate work of N.F.H. was supported by a National Science Foundation graduate fellowship and an NSF-MRL Graduate Research Assistantship for Women through the Stanford Center for Materials Research. References and Notes (1) Hansen, N. F.; Andersen, H. C. J. Chem. Phys. 1994, 101, 6032. (2) Voth, G. A.; Chandler, D.; Miller, W. H. J. Chem. Phys. 1989, 91, 7749. Voth, G. A. J. Phys. Chem. 1993, 97, 8365. (3) Garrett, B. C.; Truhlar, D. G.; Grev, R. S.; Magnuson, A. W. J. Phys. Chem. 1980, 84, 1730. Truhlar, D. G.; Gordon, M. S. Science 1990, 249, 491. (4) Messina, M.; Schenter, G. K.; Garrett, B. C. J. Chem. Phys. 1993, 98, 8525. (5) Mills, G.; Jo´nsson, H. Phys. ReV. Lett. 1994, 72, 1124. Schenter, G. K.; Mills, G.; Jo´nsson, H. J. Chem. Phys. 1994, 101, 8964. (6) McLafferty, F. J.; Pechukas, P. Chem. Phys. Lett. 1974, 27, 511. Pollak, E.; Proselkov, D. Chem. Phys. 1993, 170, 265. (7) Miller, W. H.; Schwartz, S. D.; Tromp, J. D. J. Chem. Phys. 1983, 79, 4889. (8) See, for example: Wigner, E. Trans. Faraday Soc. 1938, 34, 29. Keck, J. C. AdV. Chem. Phys. 1967, 13, 85. Anderson, J. B. J. Chem. Phys. 1973, 58, 4684. (9) Johnston, H. S. Gas Phase Reaction Rate Theory; Ronald: New York, 1966, pp 40-44. (10) Makri, N.; Miller, W. H. J. Chem. Phys. 1989, 90, 904. (11) Costley, J.; Pechukas, P. Chem. Phys. Lett. 1981, 83, 139. (12) Note that it is the exact classical rate, not the classical transition state theory rate, that is in the denominator of this expression. (13) Thirumalai, D.; Bruskin, E. J.; Berne, B. J. J. Chem. Phys. 1983, 79, 5063. (14) Jacquet, R.; Miller, W. H. J. Phys. Chem. 1985, 89, 2139.
JP951605Y
