Letter pubs.acs.org/JPCL
Semiclassical Transition-State Theory Based on Fourth-Order Vibrational Perturbation Theory: The Symmetrical Eckart Barrier John F. Stanton* Institute for Theoretical Chemistry, Department of Chemistry and Biochemistry, The University of Texas at Austin, Austin, Texas 78712, United States ABSTRACT: Semiclassical transition-state theory based on fourth-order vibrational perturbation theory (VPT4-SCTST) is applied to compute the barrier transmission coefficient for the symmetric Eckart potential. For a barrier parametrized to mimic the H2 + H exchange reaction, the results obtained are in excellent agreement with exact quantum calculations over a range of energy that extends down to roughly 1% of the barrier height, V0, where tunneling is negligible. The VPT2-SCTST treatment, which is commonly used in chemical kinetics studies, also performs quite well but already shows an error of a few percent at ca. 0.8 V0 where tunneling is still important. This suggests that VPT4SCTST could offer an improvement over VPT2-SCTST in applications studies. However, the computational effort for VPT4-SCTST treatments of molecules is excessive, and any improvement gained is unlikely to warrant the increased effort. Nevertheless, the treatment of the symmetric Eckart barrier problem here suggests a simple modification of the usual VPT2SCTST protocol that warrants further investigation.
I
and the values of Pk(E) calculated by inverting this expression to obtain θ and inserting this into eq 1 coincide with the exact quantum result.1,10,11 Together with others, Miller extended this simple harmonic SCTST approximation to second-order vibrational perturbation theory (VPT2) roughly a decade later.12,13 In VPT2,14 the independence of modes that exists in the harmonic model is lost. Hence, the theory now includes two important features: intrinsic anharmonicity of the reaction coordinate and coupling of the reaction coordinate to the orthogonal bound degrees of freedom. While not usually explicitly noted in the literature, “SCTST”, as it has been used in recent years in kinetics calculations, tends to default to this VPT2-based variant. In a number of studies, Barker and collaborators have demonstrated that this approach gives extremely good rates for simple reactions when (and only when) it is combined with potential surfaces obtained with very high-level quantum chemistry.2,3,15,16 What makes the method particularly attractive is that one need not characterize the entire potential surface associated with the reaction, or even a significant part of it. All that is required are the quadratic, cubic, and semidiagonal (iijj) quartic force constants evaluated at the transition state, which can usually be obtained for modest cost even for systems containing 5−10 atoms. Application of this theory to chemical kinetics calculations is the subject of a forthcoming review.17 While VPT2-SCTST provides excellent reaction probabilities in the energetic vicinity of the transition state, the fact that the method is based on perturbation theory means that it is not
n recent years, Miller’s semiclassical transition-state theory (SCTST)1 has proven its value in chemical kinetics studies.2−7 This theory provides a nonempirical treatment of tunneling and above-barrier reflection for levels associated with the transition state of chemical reactions, thereby providing values of Pk(E), where k identifies the particular reactive state and E is the energy. Summing Pk over all reactive states then gives the cumulative reaction probability, G(E), which is needed to calculate microcanonical and thermal rate constants.8 SCTST, and perturbative approximations to it, derive from Miller’s recognition that good action-angle variables exist at a transition state on the potential energy surface. Specifically, if the energy can be written as a Taylor series in the vibrational 1 quantum numbers nk (or, more correctly, nk + 2 ), the quantities 1
(nk + 2 ) corresponding to the bound degrees of freedom are associated with action-angle variables in the usual EBK fashion, 1 while that for the reaction coordinate (nF + 2 ) can be associated with the imaginary action iθ/π, as in WKB theory. The transmission probability, Pk(E), is then obtained from the standard semiclassical formula9 Pk(E) = [1 + exp(2θ )]−1
(1)
that is valid both below and above the barrier. For the case of a parabolic barrier with height V0 and (imaginary) harmonic angular frequency ωF, the exact expression (in atomic units) for E in terms of the classical action variables is E = V0 +
⎛
∑ ωn ⎝n⊥ + ⎜
⊥
n⊥
1 ⎞⎟ θ + iωF ⎠ 2 π © 2016 American Chemical Society
Received: June 7, 2016 Accepted: June 30, 2016 Published: June 30, 2016
(2) 2708
DOI: 10.1021/acs.jpclett.6b01239 J. Phys. Chem. Lett. 2016, 7, 2708−2713
Letter
The Journal of Physical Chemistry Letters without its shortcomings. Just as the ubiquitous VPT2 method14 for bound-state vibrational levels will perform poorly for very anharmonic potentials or for significant vibrational excitation involving motion far from the potential energy surface expansion point (a minimum), so too will VPT2-SCTST, where the expansion origin is a saddle point, misbehave for very anharmonic transition states, or in the deep tunneling regime. The former class of problematic behavior occurs already for the simple collinear H2 + H reaction;18 the latter has been emphasized by Wagner,19 who proposed a modification to improve performance at energies well below the barrier. Another rather obvious potential road to improvement is to extend the level of perturbation theory, for which this Letter provides some perspective. Because of selection rules for harmonic oscillator integrals, the VPT expansion of the energy has nonvanishing contributions only in even orders, so the next step beyond VPT2 is VPT4. While numerical calculations of bound-state vibrational energies have been carried out in the past using fourth (and higher) order VPT using numerical van Vleck perturbation theory,20,21 SCTST requires an algebraic expression for the anharmonicity constants (the so-called xij together with a constant term for VPT2; yijk 1 constants and corrections to the linear (nk + 2 ) term in the energy expression for VPT4). The latter two corrections, which are extremely tedious to evaluate, were unavailable in analytic form a few years ago and have only recently been obtained and carefully checked. Nevertheless, it is not clear that VPT4-SCTST will ever become a useful method in the field of chemical kinetics, simply because it requires all fourth, most fifth, and some sixth derivatives of the potential. These are costly to evaluate, as well as being hard to calculate with high accuracy, and the expense of a VPT4-SCTST calculation for a polyatomic molecule will be substantially greater than for VPT2-SCTST. So much so as to cause the author to be quite pessimistic about the potential of the method for seeing widespread use in future applications studies. The present paper focuses on an exploratory application of VPT4-SCTST to the energy-dependent transmission probability through a symmetric Eckart barrier.22 While just a simple onedimensional system, this potential is a venerable prototype for a reaction coordinate. Moreover, as will be seen in the following, the application of VPT4-SCTST to the Eckart barrier problem provides a solution that is in nearly perfect agreement with the exact quantum result and also suggests a simple empirical correction to much more economical (for polyatomic systems) VPT2-SCTST calculations. The latter is probably the most important aspect of this work, at least as it relates to chemical kinetics calculations. Overview of VPT and its Application to Eckart Potential. The goal of vibrational perturbation theory is to express the vibrational energy levels of a system governed by the Watson
(
Hamiltonian23 as a polynomial in the nk +
1 2
), where n
k
where the constant term G and coefficients of terms that are 1 linear, quadratic, etc. in (nk + 2 ) are expanded in a perturbation series. V0 is the value of the potential energy at the stationary point on the potential energy surface to which this expression corresponds; while usually just taken as zero and omitted from eq 3, this term has significance in the context of the present work and is retained here. It is well-known that in zeroth order, only the linear coefficients Wk (the harmonic frequencies) remain, while in second order, there is a contribution to the Xkl coefficients as well as to the constant term G, generally denoted as xkl and G0, respectively. A characteristic pattern emerges at higher orders in the expansion, as shown schematially below
where an empty square indicates that there is no contribution at that order of VPT, while a filled square indicates otherwise. That is, after a coefficient picks up its first contribution (VPT0 for Wk, VPT2 for G and Xkl, VPT4 for Yklm, VPT6 for Zklmn, etc.), additional corrections to it will occur only at orders k + 4, k + 8, and so on. Thus, in addition to providing the lowest-order Yijk terms, VPT4 also supplies a contribution to the Wk beyond the zeroth-order harmonic frequencies. Expressions for the secondorder contributions G0 and xkl are well-known and can be found in the literature;14,24−26 checking of correct general (multimode) algebraic expressions for the fourth-order contributions to Wk and Yijk is nearly complete, and these will be published elsewhere. However, for the simple case of a one-dimensional oscillator, the relevant equations are given below, where the energy, as is customary in VPT, is given in wavenumber units (E/hc) and ν̃e is the harmonic wavenumber: W = νẽ −
Y=
is the
+
⎛
∑ Yklm⎝nk + ⎜
klm
288
9216νẽ
−
+
7ϕ3ϕ5 288νẽ
5ϕ6 1152
−
17ϕ4 2 2304νẽ
(4)
+
25ϕ32ϕ4 384νẽ 2
−
235ϕ3 4 6912νẽ 3 (5)
quantum number of the kth mode in the zeroth-order (harmonic oscillator) picture, viz. ⎛ 1⎞ E = V0 + G + ∑ Wk ⎜nk + ⎟ + ⎝ 2⎠ k
ϕ6
67ϕ4 2
while the second-order parts of G and X are given by
⎛ 1⎞ ∑ Xkl⎛⎜⎝nk + 1 ⎞⎟⎜ nl + ⎟ ⎠ ⎝ 2 2⎠ kl
X=
1 ⎞⎟⎜⎛ 1 ⎞⎛ 1⎞ nl + ⎟⎜nm + ⎟ 2 ⎠⎝ 2 ⎠⎝ 2⎠
G=
⎛ 1 ⎞⎛ 1 ⎞⎛ 1 ⎞⎛ 1⎞ + ∑ Zklmn⎜nk + ⎟⎜nl + ⎟⎜nm + ⎟⎜nn + ⎟ + ··· ⎝ ⎠⎝ ⎠⎝ ⎠⎝ 2 2 2 2⎠ klmn
ϕ4 16 ϕ4 64
−
−
5ϕ32 48νẽ
(6)
7ϕ32 576νẽ
(7)
In the equations above, the ϕn are nth derivatives of the energy in wavenumbers with respect to the dimensionless normal coordinates q, and of course ν̃e = ϕ2. For a one-dimensional
(3) 2709
DOI: 10.1021/acs.jpclett.6b01239 J. Phys. Chem. Lett. 2016, 7, 2708−2713
Letter
The Journal of Physical Chemistry Letters potential governing a particle of mass m, the ϕn are related to geometrical derivatives of the potential f n [ fn ≡ ϕn =
f
dnV ] dx n e
( )
n n n /2 n /2 γ m νẽ hc
fourth-order VPT for the Eckart potential, due to a cancellation similar to that beyond VPT2 for all of the constants of eq 3 for the Morse oscillator, where the simple second-order theory is exact. VPT and SCTST. In semiclassical transition-state theory, the probability that a reactive state k at energy E proceeds to products is given by the usual expression
via
(8)
where
⎛ c ⎞1/2 γ ≡ 2π ⎜ ⎟ ⎝h⎠
Pk(E) = [1 + e 2θk(E)]−1 (9)
θk(E), which is analogous to the barrier penetration integral, is determined from the polynomial energy expression, eq 3, as applied to a transition state rather than to a potential energy surface minimum. In such a case, one of the normal modes (designated as F) corresponds to the reaction coordinate and has 1 an imaginary harmonic frequency; the corresponding nF + 2
In the context of SCTST, another consideration becomes important. Specifically, when the harmonic wavenumber ν̃e is imaginary, it is convenient to define the corresponding dimensionless normal coordinate in the usual way, but based on the
(
magnitude |ν̃e|; one then attaches a phase factor exp −
ikπ 4
(
) to the
⎛ 1 ⎞⎟ iθ ⎜n + ⇒ ⎝ F 2⎠ π
≈ 29.347 (where α is the fine structure
2π α
f n(a.u.) n + 2/2
( 2απ )
mn /2νẽ n /2
(10)
For the symmetric Eckart potential governing the motion of a particle with mass m V = V0 sech2(x /a)
G (E ) =
(11)
f4 =
(12)
16V0 a4
f6 = −
(13)
E = V0 + |νẽ |i
272V0 a6
(14)
so that
where the odd derivatives vanish because the potential is symmetric. It is easily shown that the coefficients of the VPT4 polynomial given by eq 3 (with E in atomic (bohr−1) wavenumber units) become (with ℏ = 1)
θ=
G=
2 −1 −1 ⎡ ω b ⎤ ⎥ ⎢ = 2 2πc ⎣ 16V0 ⎦ 16ma cπ
⎡ 1/2 1 V W= ⎢ 0 − cπ ⎢⎣ 2ma2
⎤ 3 ⎤ ⎡ ⎥i = 1 ⎢ω − ω b ⎥i b 2πc ⎣ 32V0 2 ⎦ 512V0m3a6 ⎥⎦
(18)
(23)
θ θ2 −X 2 π π
(24)
which has the solutions
(16)
Y=0
(22)
π (V0 − E) |νẽ |
E = V0 + G − |W |
1
(17)
|ν ̃ | iθ = V0 − e θ π π
which, when inserted into eq 19, gives the exact quantum result for P(E) associated with a parabolic barrier.10 For VPT2, the corresponding energy equation is
(15)
2 −1 −1 ⎡ ω b ⎤ ⎥ = 4G0 ⎢ X= = 2 2πc ⎣ 4V0 ⎦ 4ma cπ
(21)
The description here gives a brief overview of how VPT-SCTST is used in practice. VPT-SCTST for the Symmetric Eckart Barrier. In zeroth-order (VPT0), the prescription above applied to a one-dimensional oscillator gives
2V0 a2
∑ Pk(E) k
the geometric derivatives ultimately needed for a VPT4 treatment are f2 = −
(20)
As made implicit earlier, when inserted into eq 3, this provides an equation that can be inverted to provide θ(E) for a reactive state at energy E, where the state is described by quantum numbers nk corresponding to the modes orthogonal to the reaction coordinate (k ≠ F). The cumulative reaction probability G(E), fundamental to the calculation of both microcanonical (constant E) and thermal rate constants is then given by summing eq 19 over all reactive quantum states
constant, c = α−1 in atomic units), and eq 8 can be written as ϕn(a.u.) =
)
terms in the energy expression are written (by analogy with the WKB method) in terms of the imaginary actions
force constants associated with the kth derivative of the potential with respect to the unstable dimensionless coordinate. For a potential given in atomic units, with mass and harmonic wavenumber ν̃e also in atomic units (me and bohr−1, respectively), γ =
(19)
θ±(E) = π
−|W | ∓ [|W |2 + 4(V0 + G − E)X ]1/2 2X
(25)
where θ+(E) is chosen because it connects with the VPT0 result in the limit that the second-order quantities G and X vanish. For the general VPT4-SCTST case, the equation for θ is cubic, but the coefficient of the cubic term (Y) vanishes for the Eckart potential. Hence, eq 25 is again operative here, with the proviso that W now contains the fourth-order contribution in eq 16. Inserting eqs 15−17 into those above gives the following results
where the magnitude of the imaginary barrier (angular) frequency is denoted here and in the following by ωb. It should be noted that the cubic coefficient Y vanishes identically in 2710
DOI: 10.1021/acs.jpclett.6b01239 J. Phys. Chem. Lett. 2016, 7, 2708−2713
Letter
The Journal of Physical Chemistry Letters in terms of the barrier height, V0, and angular frequency, ωb, for the Eckart barrier problem θ VPT0‐SCTST =
θ VPT2‐SCTST
πV0 ⎧ E⎫ ⎨1 − ⎬ ωb ⎩ V0 ⎭
1/2 ⎫ ⎧ ⎡⎛ E ⎞ 2πV0 ⎪ ω b2 ⎤ ⎪ ⎥ ⎬ ⎨1 − ⎢⎜ ⎟ + = ⎢⎣⎝ V0 ⎠ 16V0 2 ⎥⎦ ⎪ ωb ⎪ ⎩ ⎭
θ VPT4‐SCTST =
(26)
(27)
1/2 ⎫ ⎧ ⎡⎛ E ⎞ ω b2 ω b4 ⎤ ⎪ 2πV0 ⎪ ⎢ ⎥ ⎨1 − ⎬ − + ⎜ ⎟ ⎢⎣⎝ V0 ⎠ 1024V0 4 ⎥⎦ ⎪ ωb ⎪ 32V0 2 ⎩ ⎭
(28)
For the Eckart potential studied by Miller and Seideman (m = 1061, V0 = 0.0156, a = 0.734, all in au,29 for which ωb = 0.00738794 au), which provides a model for the collinear H + H2
Figure 2. Similar to Figure 1, but displaying the deep tunneling regime.
is given quite accurately by VPT2-SCTST (note that this shift in energy is given here by the value of G). VPT0-SCTST, however, while appropriate for a parabolic barrier, is decidedly worse than VPT2-SCTST and WKB, but still gives a plausible treatment of tunneling near the barrier, which is the most important region for most chemical kinetics studies. It is, however, completely useless in the deep tunneling regime, as one should expect. What is truly remarkable here is the quality of the VPT4SCTST results. At 0.7V0, VPT4-SCTST and exact Eckart transmission probabilities differ by one part in 104, and the agreement falls to only one part in 103 at 0.2V0, where P(E) < 10−6. Indeed, one has to go so deep into the tunneling regime (where P(E) ≈ 10−10) to find a 1% error for VPT4 that any implication for associated chemical kinetics calculations is negligible. The excellent performance of VPT4-SCTST for the Eckart barrier motivates investigation of specifically how it deviates from the exact result. The starting point for this analysis is to recognize from eqs 19 and 29 that the SCTST and exact tunneling probabilities for VEckart will become equal if
Figure 1. Ratio of barrier transmission coefficients computed by WKB and VPT-SCTST approaches to the corresponding exact quantum result for the Eckart barrier described in the text, for energies in the vicinity of the bare adiabatic barrier. The blue line is VPT0-SCTST (the parabolic barrier approximation); pink is WKB, red VPT2-SCTST, and black VPT4-SCTST.
exp(2θ ) =
reaction, results for P(E) are shown in Figures 1 and 2, where the exact quantum result27 1 P(E) = ⎛ 2 2 2 ⎞ 1+
cosh2⎜ ⎝
4π (16V0 − ω b ) 16ω b2
⎛ sinh2⎜ ⎝
4π 2V0E ω b2
⎞ ⎟ ⎠
⎛ cosh2⎜ ⎝
4π 2(16V0 2 − ω b2)
⎛ sinh2⎜ ⎝
16ω b2 4π 2V0E ω b2
⎞ ⎟ ⎠
⎞ ⎟ ⎠
(30)
or ⎧ 2 2 2 ⎞⎫ ⎛ cosh⎜ 4π (16V0 −2 ω b ) ⎟ ⎪ ⎪ ⎪ 16 ω ⎝ ⎠⎪ b ⎬ θ = ln⎨ ⎛ 4π 2V0E ⎞ ⎪ ⎪ ⎟ sinh⎜ ⎪ ⎪ ⎝ ω b2 ⎠ ⎭ ⎩
⎟ ⎠
(29)
is compared to the three SCTST methods described earlier, as well as to the corresponding WKB calculation (as generalized by the uniform semiclassical approximation to above-barrier energies). It should be noted that the performance of VPT2-SCTST is excellent, as has been established elsewhere.18 Near and above the barrier, VPT2-SCTST is close to the exact quantum result (
