DOI:
10.1039/C7CS00602K
(Review Article)
Chem. Soc. Rev., 2017,
46, 75487596
Variational transition state theory: theoretical framework and recent developments
Received
14th August 2017
First published on 22nd November 2017
This article reviews the fundamentals of variational transition state theory (VTST), its recent theoretical development, and some modern applications. The theoretical methods reviewed here include multidimensional quantum mechanical tunneling, multistructural VTST (MSVTST), multipath VTST (MPVTST), both reactionpath VTST (RPVTST) and variable reaction coordinate VTST (VRCVTST), systemspecific quantum Rice–Ramsperger–Kassel theory (SSQRRK) for predicting pressuredependent rate constants, and VTST in the solid phase, liquid phase, and enzymes. We also provide some perspectives regarding the general applicability of VTST.
Junwei Lucas Bao
 Junwei Lucas Bao joined the Department of Chemistry at the University of Minnesota Twin Cities in 2013, and he is currently pursuing a PhD degree in theoretical chemistry under the supervision of Professor Donald G. Truhlar. Currently, his main research interests include: (1) multistructural and multipath variational transition state theory with multidimensional tunneling, theory of pressuredependent rate constants, torsional anharmonicity, and their applications in combustion and atmospheric chemistry; (2) multiconfiguration pairdensity functional theory (MCPDFT), including its separatedpair (SP) version with the correlatingparticipatingorbital (CPO) scheme, and their applications in modeling transition metal compounds and spinsplittings of divalent radicals. 
Donald G. Truhlar
 Don Truhlar is native of Chicago. He received a PhD from Caltech and has been on the Chemistry faculty at the University of Minnesota since 1969, now as Regents Professor. He was honored to receive the ACS Award for Computers in Chemical and Pharmaceutical Research, NAS Award for Scientific Reviewing, ACS Peter Debye Award for Physical Chemistry, WATOC Schrödinger Medal, Dudley Herschbach Award for Collision Dynamics, RSC Chemical Dynamics Award, and APS Earle Plyler Prize for Molecular Spectroscopy and Dynamics. His current research interests center on dynamics and electronic structure, especially photochemistry, transition metal catalysis, and density functional theory. 
1. Introduction
Chemical kinetics has a long history in chemical science. In 1864, Waage and Guldberg formulated the law of mass action,^{1} which can be considered as the first published work developing the theory of chemical kinetics. About 25 years later, Arrhenius proposed his well known empirical relation for the temperature dependence of reaction rates.^{2} Both of these studies still exert a profound influence on modern chemical kinetics. In the 1920s, the concept of a transition state emerged, and in the 1930s it became a central theoretical foundation for quantitative theory, especially because of the work of Eyring, Wigner, Polanyi, and Evans.^{3} The direct observation of transition state structures by various spectroscopic techniques is considered to be one of the major challenges in modern chemical science.^{4}
For a long time, transition state theory was viewed somewhat skeptically as far as its relevance to quantitative work. For example, a group of distinguished researchers wrote in 1983 that “The overall picture is that the validity of the transition state theory has not yet been really proved and its success seems to be mysterious.”^{5} This began to change though when accurate quantum mechanical calculations became available for prototype reactions of small molecules, and when accurate semiclassical tunneling methods became available.^{6} Now transition state theory is understood on a firm basis, especially since the observation of quantized transition state level structure in converged quantum mechanical dynamics calculations^{7–9} has provided a solid quantum mechanical foundation for transition state theory in terms of shortlived resonance states.^{10–13}
Transition state theory is basically a theory for electronically adiabatic reactions, i.e., reactions that occur on a single Born–Oppenheimer potential energy surface. There are some extensions to multiplesurface reactions, but they involve additional assumptions that are not traditionally considered to be part of transition state theory, and we will not consider them here. Furthermore, transition state theory in its basic form predicts reaction rates in equilibrated ensembles, in particular canonical reaction rates in thermal ensembles and microcanonical reaction rates in ensembles with fixed total energy or fixed total energy and total angular momentum; thus it is not designed for statespecific reaction rates or for predicting product internalstate distributions, although sometimes it is applied to such properties by making additional assumptions.
Transition state theory is most powerful in the form of variational transition state theory^{14–25} (VTST), especially when combined with modern electronic structure theory, multidimensional tunneling methods, and statistical mechanical and quantum mechanical theory for anharmonicity, and it has become a powerful tool in understanding and predicting chemical reaction mechanisms and reaction rates.^{26–36}
We have reviewed VTST elsewhere in a variety of ways,^{6,20–25,30,33,35,37–41} but the current review has a special focus on several modern enhancements that were not available when those earlier reviews were written. In Section 2, we cover the fundamental principles of variational transition state theory in the context of gasphase reactions; this section uses a reaction coordinate associated with a reaction path, and it includes tunneling. The next four sections continue the treatment of gasphase reactions, covering multistructural transition states (Section 3), variable reaction coordinates (Section 4), pressure effects (Section 5), and tests of the applicability of the theory (Section 6).
Sections 1–3 apply to tight transition states, while Section 4 applies for loose transition states. Section 7 discusses extension of the treatment to condensedphase reactions, and Section 8 discusses selected applications. We conclude this review in Section 9.
Most of the computational methods reviewed here are available in the Polyrate computer code.^{42}
2. Fundamentals of gasphase variational transition state theory
A detailed pedagogical review of the fundamental theory underlying gasphase variational transition state theory, including derivations, is available elsewhere.^{41} In this section we present a briefer summary including a few new developments. We start by using classical mechanics to derive variational transition state theory for gasphase reactions (Section 2.1), and then we introduce quantum effects into VTST, in particular vibrational mode quantization in the partition functions (Section 2.2) and quantum mechanical tunneling in the treatment of reactioncoordinate motion, including multidimensional coupling of reaction coordinate motion to other coordinates (Section 2.3). We emphasize the importance of utilizing curvilinear coordinate for carrying out the VTST calculations in Section 2.5. Stateselected extensions of VTST are presented in Section 2.6.
2.1 Classical variational transition state theory
Transition state theory is first justified in classical mechanics; then we add quantization effects. In classical mechanics, we work with trajectories. At any given time, a trajectory is a point in phase space, which is the space of all the coordinates and their conjugate momenta. A point in phase space is called a phase point.
Transition state theory may be applied to both thermal ensembles and microcanonical ensembles. In either case the first assumption is that the reactants are in local equilibrium. For thermal ensembles, this means that the population of reactant states is given by a Boltzmann distribution, even if reactants are not in equilibrium with products. For a microcanonical ensemble, this means for classical systems that all cells of equal volume in reactant phase space are equally populated; for microcanonical quantum systems it means that all reactant states of a given energy are equally populated. In the rest of this subsection we assume a classical mechanical world.
The maintenance of reactant equilibrium means that the rate of thermalizing collisions for maintaining the equilibrium states of these species is at least as fast as the rate at which these states are depopulated by the chemical reaction. Since the phase points in the reactant region are assumed to be in thermal equilibrium, Liouville's theorem guarantees that as the time evolves, the Boltzmann distribution is also satisfied in the transition state region except that if products are absent then phase points corresponding to trajectories that originated in the product region are missing.^{43,44} This localequilibrium assumption in the reactant region is assumed valid for most gasphase bimolecular reactions and reactions in the liquid phase because the inelastic collisions of the reactants are efficient enough to repopulate their rovibrational states and maintain local equilibrium in phase space, but it is not applicable for unimolecular reactions (and some bimolecular reactions) in the gas phase at low pressure. When the equilibrium assumption breaks down, “thermal” rate constants become pressure dependent, and the thermal rate constants calculated by using the transition state theory are the highpressure limit rate constants. In the highpressure limit, the word “thermal” refers to rate constants averaged over Boltzmann distributions of both reactants and bath gases; but at lower pressures, for reactions whose rate constants are pressuredependent, “thermal” refers to a Boltzmann distribution of bath gases, but there is a pressure and temperaturedependent nonBoltzmann energy distribution of reactants. Thus, in the falloff region, rate constants are functions of both pressure and temperature. We will discuss pressure dependence in Section 5.
It is convenient to introduce massscaled coordinates,^{22,45} which may also be called isoinertial coordinates. With such a coordinate system, the motion of an Natom system in three dimensions can be reduced to that of a singleparticle with a reduced mass μ moving in 3Ndimensional space. In other words the reduced mass becomes the same for all directions of motion (and hence the word “isoinertial”). The reduced mass μ is usually chosen to be 1 atomic mass unit, or it can be set to be the reduced mass of the relative translational motion for a bimolecular reaction (one can use any value as long as one is consistent). The value of μ sets the scale for the reaction coordinate s; however, the value of any physical observable is independent of the value of μ. If one sets μ = 1 amu, the numerical value of s in Å is the same as the numerical value of s in amu^{1/2} if one used the familiar massweighted coordinates^{46} of spectroscopists. The motions of all the particles are governed by a potential energy surface V({r_{i}}) that depends on the set of N threedimensional positions r_{i} of the particles, each with its own mass. Equivalently speaking, the motion may be described as that of a single particle with reduced mass μ governed by V(q), where q is a 3Ndimensional generalized massscaled coordinate, and its conjugate momentum is p.
Now we present classical transition state theory. Consider an ensemble of systems in phase space. Each phase point (q,p) in the 6Ndimensional phase space represents an Natom system, and the motion of each point is governed by the following set of Hamilton equations:

 (1) 
where the classical Hamiltonian is

 (2) 
Liouville's theorem states that the phase space volume is conserved as time evolves; it can be written in the form of a continuity equation as follows:

 (3) 
where
∇ is the 6
Ndimensional divergence operator,
v is a 6
Ndimensional velocity of a phase point, and
ρ is the density of points in phase space. Integrating
ρ over the whole phase space yields the total number of classical phase points (
i.e., the total number of systems in the ensemble).
The number of systems located in the reactant region is N^{R}, and it can be evaluated by integrating the density of phase points over the phase space of the reactant:

 (4) 
where
R indicates that the volume integral is carried out over the reactant region. Differentiating
eqn (4) with respect to
t and plugging into
eqn (3) yields

 (5) 
Then we can use Gauss's theorem to rewrite the volume integral in
eqn (5) as a surface integral,

 (6) 
where
n is a unit vector normal to the surface
S pointing out of the volume
R, and
S is a (6
N − 1)dimensional hypersurface that separates the reactant region from the product region; this hypersurface is called the dividing surface.
The local oneway flux of systems going from reactants to products through this dividing surface is

 (7) 
where we consider the portion,
S^{+}, of the dividing surface for which the flux is positive. We define the coordinate
q_{3N} as the coordinate normal to surface
S, and we call this coordinate the local reaction coordinate
z. Note that the surface
S is locally normal to some global reaction path from reactants to products. The dividing surface is put at the location where the global reaction coordinate is
s*, and the local reaction coordinate is
z*. The local oneway flux is evaluated as:

 (8) 
As we mentioned at the beginning of this subsection, TST assumes the reactants are in local equilibrium at temperature T and the density of classical phase points satisfies the Boltzmann distribution
where
H is the Hamiltonian,
k_{B} is Boltzmann constant,
T is temperature, and
ρ_{0} is a constant. Plugging
eqn (9) into
eqn (8) and (4) yields

 (10) 

 (11) 
Now, we define a generalized transition state (GT) as a (6N − 2)dimensional fixedz hyperplane with the following Hamiltonian H^{GT}:

 (12) 
The generalizedtransitionstate Hamiltonian
H^{GT} is defined by removing the kinetic energy associated with the motion of the local reaction coordinate
z from the total Hamiltonian, and fixing the value of
z at some specific value
z*. Notice that
z is a parameter rather than a variable in
H^{GT},
i.e.,
H^{GT} depends parametrically on
z. Substituting
eqn (12) into
eqn (10), we obtain the oneway flux
F^{GT}(
T,
z*) through a generalized transition state fixed at
z = z*, as

 (13) 
where GT indicates that the integral is carried out over the generalized transition state.
The local equilibrium oneway flux F^{GT}(T,z*) gives an upper bound to the global reactive flux F(T) if we arbitrarily put our dividing surface at the saddle point in this classical dynamics approach. There is no guarantee that the trajectories in phase space that are originating from the reactant region toward the product region will cross the dividing surface once and only once. In the conventional transition state theory, one simply puts the dividing surface at the saddle point and makes a nonrecrossing assumption that there is no recrossing on this surface.
In order to have the best estimation of the classical rate constant, the ideal location to put the dividing surface S would be the location at which F^{GT}(T,z*) = F(T), which means all the trajectories originating from the reactant region cross the surface only once (no recrossing) and end up in the product region. If one is willing to make the transition state hypersurface arbitrarily convoluted in shape, or if one considers simple enough systems at low energy,^{47} one can make classical transition state theory exact. For real potential energy surfaces and dividing surfaces that are simple enough to allow for the calculation of the flux, completely eliminating all the recrossing is almost impossible, so we would put the dividing surface at the socalled dynamical bottleneck, which is the position at which the recrossing is minimized. Classical transition state theory will overestimate the classical rate constant, but the variational procedure will minimize this overestimation.
For a bimolecular reaction (A + B), the classical generalized transition state theory rate constant k^{GT}_{Classical}(T), which is the flux per volume divided by the product of the concentration of reactants, obtained by the above procedure is

 (14) 
where
V is the threedimensional realspace volume of the system, and GT indicates that the integral is carried out over the generalized transition state. We can use the following definition of the classical partition function for
A 
 (15) 
where
h is Planck's constant, with an analogous expression for B. Recall that Planck's constant occurs in the classical partition function because we require the classical partition function to equal the quantum one in the classical limit.
^{48} We also define a classical partition function per volume as

Φ_{Classical} = Q_{Classical}(T)/V  (16) 
For the generalized transition state, we freeze one degree of freedom (the reaction coordinate) and take the local zero of energy as the classical potential energy V_{RP} evaluated at z = z* on the reaction path (RP). The generalized transition state classical partition function per volume for the transition state becomes

 (17) 
The total Hamiltonian can be written as the sum of the Hamiltonians of the reactants,
H_{A} and
H_{B}, and using

N^{R}_{X} = ρ^{X}_{0}h^{3NX}VΦ^{X}_{Classical}  (18) 
where
N^{R}_{X} is the number of atoms in X (
N^{R}_{A} +
N^{R}_{B} =
N), and the constant
ρ_{0} =
ρ^{A}_{0}ρ^{B}_{0}, we obtain

 (19) 
in which
X is A or B,
X represents that the volume integral is carried out over reactant region, and
H^{X} is the Hamiltonian of the reactant X.
We then use eqn (15)–(18) to rewrite eqn (14) as

 (20) 
Follow a similar derivation, the classical generalized TST rate constant for unimolecular reactions is:

 (21) 
A convenient choice for the reaction path is the minimum energy path (MEP) in isoinertial coordinates.^{49} The global reaction coordinate s is then the signed distance along the MEP. An MEP is the union of the paths of steepest descent in isoinertial coordinates from the saddle point towards reactants and products. The physical meaning of the MEP is that it is the path followed by a classical trajectory that starts at the saddle point but is continuously damped so as to always have only an infinitesimal velocity. A positive value of s represents for the direction toward the product; s = 0 is the saddle point; and a negative value of s represents the direction toward the reactant side. Starting from the saddle point with coordinates x^{‡} (which may also be called x(s = 0)), the geometry of a nonstationary point that is close to the saddle point on the reaction path is computed using:

x(s_{1} = ±δs) = x(s = 0) ± δs·L(x^{‡})  (22) 
where δ
s is the step size (typically chosen to be 0.001 to 0.003 Å) along the MEP, and
L(
x^{‡}) is the normalmode eigenvector associated with the imaginary frequency (it is a column of the orthogonal matrix that diagonalizes the Hessian matrix evaluated at the saddle point). If we start from any nonstationary point
x(
s_{j}) on MEP, then the next point can be computed by using the normalized gradient as:

 (23) 
At each value of s, we can define a new coordinate system by rotating and translating the original coordinates so that the coordinate z that we mentioned in the generalized TST is tangent to MEP at s, and any surface that is orthogonal to MEP can be described by z* = 0. The dividing surface intersects MEP at the reaction coordinate s*. Note that s* depends on temperature.
In canonical variational transition state theory (CVT), we variationally optimize the position of the dividing surface to minimize the generalized transition state rate constant k^{GT}_{Classical}(T) since the local oneway flux provides an upper bound of the global reactive flux. The CVT rate constant is computed by

 (24) 
Therefore

 (25) 
which leads to the following CVT rate constant for bimolecular reactions

 (26) 
and for unimolecular reactions

 (27) 
Another way to understand the canonical variational approach is to rewrite the classical generalized TST rate constant in terms of the standardstate Gibbs free energy of activation Δ
G^{GT,0}_{act}(
T,
s).

 (28) 
where
K^{0} is the reaction quotient (for the formation of TS) evaluated at the standardstate concentration of each species:

 (29) 
where, for gasphase bimolecular reaction with unit cm
^{3} molecule
^{−1} s
^{−1}, (
c°)
^{−1} is
k_{B}T/
p°, where
p° is 1 bar; for liquidphase reaction, the standardstate concentration is 1 mole L
^{−1}.
CVT is equivalent to maximizing the free energy of activation

 (30) 
In a microcanonical ensemble (NVE ensemble), the microcanonical energyresolved classical rate constant (which is obtained by firstly integrating the E, Jresolved rate constant k^{GT}(E,J,s) over the angular momentum J distribution) is

 (31) 
where
N^{GT}(
E,
s) is the phasespace volume of the generalized transition state (it is, quantum mechanically, the number of states for the generalized transition state), and
ρ_{R}(
E) is the density of states for reactant. The density of states is the inverse Laplace transform of the partition function. In microcanonical variational transition state theory (μVT), the variational rate constant is obtained by minimizing
eqn (31) with respect to
s:

 (32) 
and therefore the location of microcanonical variational transition state is energydependent. To obtain canonical μVT rate constant, one needs to perform the following integration:

 (33) 
where
P_{R}(
E) is the normalized population distribution of the reactant:

 (34) 
where
Q_{R}(
T) is the canonical partition function of reactant. If the dividing surface is located at the saddle point of the potential energy surface, then
eqn (33) is equivalent to the rate constant expression of the conventional transition state theory. It can be shown that

k^{CVT}(T) ≥ k^{μVT}(T)  (35) 
Inserting eqn (34) into (33) shows that the canonical rate constant may be considered as a Laplace transform of the microcanonical one, or the microcanonical rate constant may be considered to be an inverse Laplace transform of the canonical one. This will be a useful perspective in Section 5.
2.2 Quantization: quasiclassical variational transition state theory
The classical canonical variational transition state theory cannot be rigorously quantized because in order to define the generalized transition state appropriately, we need to know the local reaction coordinate z and its conjugate momentum p_{z} simultaneously, which is a violation of Heisenberg's uncertainty principle. However, the quantum mechanical effects along the reaction path usually cannot be ignored, and therefore we put in the quantum effects approximately to obtain what may be called quantized variational transition state theory. In order to quantize classical variational transition state theory,^{38,50–56,77} we need to: (1) quantize the translational, rotational, and vibrational degrees of freedom (the electronic degrees of freedom have already been quantized by the Born–Oppenheimer approximation); (2) quantize the corresponding partition functions; and (3) quantize the motion along the reaction coordinate, which includes quantum mechanical tunneling (including nonclassical reflection). We will discuss quantum mechanical tunneling in Section 2.3.
First, we quantize all the degrees of freedom except for the vibrational mode associated with the reaction coordinate s at a given generalized transition state. This degree of quantization was called hybrid in early papers on VTST; now we call it quasiclassical. In many cases we ignore the coupling between some of the degrees of the freedom; that is, we assume independent normal modes of vibration, and we neglect vibration–rotation coupling. VTST theory is not restricted to neglecting these couplings, but computations of the partition functions are more complicated if they are included. (We do include such coupling in some of the approximations considered later in this section and in Sections 3.1 and 4.1.)
The quantized translational motion is described by a particle in 3dimensional box; since box length is macroscopic, this agrees exactly with the classical result. For overall rotation, usually, it is sufficient to use the classical rotor approximation, provided one also correctly includes the rotational symmetry number (which is a quantum effect); nevertheless, at accessible temperatures including 300 K and below, all of the diatomic hydrides, many important free radicals (e.g., OH, CH, and CH_{2}), and a few polyatomics have large rotational constants, for which the classical approximation is not always sufficiently accurate. For H_{2}, it is advisable to always use the quantal rotational partition function.
For vibrational motion, the traditional procedure is to use the harmonic oscillator approximation, which assumes a quadratic potential. A quadratic potential has only one minimumenergy structure, whereas real molecules often have multiple structures (conformations) due to torsions and/or rings. In some cases, the singlestructure approximation causes large errors, and eliminating these errors is the subject of Section 3.
If the goal is a quantitative calculation, one should explicitly consider anharmonicity even for stretches and bends in a single structure (this vibrational anharmonicity needs to be corrected because the harmonicoscillator approximation truncates the series expansion of the potential energy at the secondorder terms). Vibrational anharmonicity can be included by vibrational selfconsistentfiled theory (VSCF), vibrational perturbation theory (VPT), vibrational configurational interaction (VCI) and other VSCFbased methods,^{57} but this can get expensive. To avoid the cost of those higherlevel vibrational calculations, it is convenient use a scale factor^{58} (for a given level of electronic structure theory) to scale all the computed harmonic vibrational frequencies in order to account for vibrational anharmonicity; such scale factors can also partially correct for systematic errors in the chosen level of electronic structure theory. Using the harmonic oscillator formulas with effective frequencies (such as those obtained by scaling) is called the quasiharmonic oscillator (QHO) approximation. (Sometimes it is just called the harmonic oscillator approximation, but that can lead to confusion.) Note that the scaled frequencies can (and usually do) account for vibrational mode coupling as well as intramode anharmonicity, so the QHO approximation is a way to include nonseparability in a separablemode formalism.
For the purpose of thermochemical kinetics computations, the scale factor we often choose is the one that is developed for reproducing an accurate zero point energy; this is denoted as λ^{ZPE}, and it includes the corrections both for systematic errors in the electronic structure level and for the vibrational anharmonicity. The scale factors^{58} are developed based on a database that consists of 15 representative small molecules, whose experimental harmonic vibrational frequencies, fundamental frequencies and zeropoint vibrational energies are available. The ZPE scale factor (λ^{ZPE}) can be written as the product of the harmonic frequency scale factor (λ^{HO}, which corrects the systematic error of a given electronic structure method) with the vibrational anharmonicity scale factor (λ^{anh}, which corrects the vibrational anharmonicity).
Instead of using the generic scale factors, the vibrational anharmonicity scale factor for a given structure can be computed by using hybrid,^{59} degeneracycorrected,^{60} secondorder^{61–63} vibrational perturbation theory (HDCVPT2); and it is equal to the ratio of HDCVPT2 computed anharmonic ZPE to the harmonicoscillator ZPE. (The use of the HDCVPT2 method, like the use of vibrational scale factors, includes the effects of mode–mode coupling.)
To properly take account of anharmonicity, we have usually assumed that it is sufficient to include anharmonicity in all bends and stretches by a single scaling factor. Note that the scale factors designed to treat the zeropoint energy are primarily dictated by the stretching frequencies since they are larger than the bending and torsion frequencies and hence their contributions dominate the zeropoint energy. Since stretches usually have negative anharmonicity, frequencyscaling factors are usually smaller than unity, which means that partition functions are increased by using a frequencyscaling factor. However, anharmonicity of a bend mode could either increase or decrease the partition function. To achieve the most accurate possible thermodynamic and kinetic calculations requires treating anharmonicity in each highfrequency mode/nontorsional mode individually. One approximation for computing such vibrational partition function with vibrational anharmonicity is called simple perturbation theory (SPT);^{64} it only requires calculation of the fundamental frequencies, which can be computed by the abovementioned HDCVPT2 method.
The quantized total partition function of species X can be then written as

Q(X) = Q_{trans}(X)Q_{rot}(X)Q_{vib}(X)Q_{elec}(X)  (36) 
with the vibrational partition function in the QHO approximation given by

 (37) 
where for the reactant and product, the vibrational degree of freedom
F = 3
N − 5 for linear molecule, or 3
N − 6 for nonlinear molecule, 3
N − 6 for linear transition states, and 3
N − 7 for nonlinear transition states. And
ν_{i} is the vibrational frequency for the normal mode
i.
Note that we have introduced the zero point energy, we have to specify what we take as the zero of energy for partition functions. Textbooks routinely take the zero point level as the zero of energy; then partition functions tend to unity (or, more generally, to the degeneracy of the ground state) when the temperature goes to zero. In our work (and in this whole article) the zero of energy for the vibrational partition function is chosen to be V_{MEP}(s). We also define the zero of energy for the electronic partition function as V_{MEP}(s), and hence, without considering spin–orbit coupling or lowlying electronic states, the electronic partition function is just equal to the degeneracy of the ground electronic state, but the vibrational partition function is not (and we do include lowlying electronic states in the electronic partition function when they are present – see the next paragraph). The rotational partition function (which is approximated by the classical rotor approximation) is determined by the principal moments of inertial of the species. For the generalized transition state, Q^{GT} is a function of reaction coordinate s.
In thermochemical kinetics we usually assume the molecules are in their ground electronic states. The zero of energy is set equal to be electronic energy of the ground state, the electronic partition function is

Q_{elec} = g_{0} + g_{1}e^{−E1/kBT} +⋯  (38) 
where
g_{0} and
g_{1} are the degeneracy of ground state and first excited state, and
E_{1} is the adiabatic excitation energy of the first excited state (
i.e., the first excited state energy relative to the ground state). Usually the first excited state is much higher than the ground state, and thus only including the first term is usually sufficient,
i.e., the electronic partition function is equal to the groundstate degeneracy. Notable exceptions are reactions involving halogen atoms and the OH radical.
Next we need to consider the vibrational zeropoint energy (ZPE), which is of central importance in the interpretation of the kinetic isotope effects.^{65–74} The generalized normal mode analysis is carried out along the MEP, and the total zeropoint energy at each s is computed as the summation of the ZPEs of all the vibrational modes that are orthogonal to the reaction coordinate. We define the vibrationally adiabatic groundstate potential energy V^{G}_{a} as follows:

V^{G}_{a} = V_{MEP}(s) + ε^{G}(s)  (39) 
where
ε^{G}(
s) is the zeropoint energy of the transverse vibrational modes. The generalized normal mode analysis along the reaction coordinate is carried out with curvilinear coordinates rather than with rectilinear coordinates or Cartesian coordinates. In order to eliminate or minimize the appearance of unphysical imaginary frequencies for modes transverse to the reaction coordinate, one should use curvilinear coordinates, which are physically more meaningful.
^{75,76} Note that only at the stationary points are the harmonic vibrational frequencies independent of the choice of coordinates.
One important practical issue about the formulas for computing rate constants is the choice of the zero of energy for the calculation, which must be done consistently. If one chooses the zero of energy as the bottom of the potential energy well, i.e., the equilibrium potential energy, for the lowestenergy structure of the reactant in a unimolecular reaction, or as the summation of the equilibrium potential energies of the two (infinitely separated) reactants in a bimolecular reaction, then the vibrational partition function must be computed accordingly, and this is why we include e^{−hνi/2kBT} in the numerator of eqn (37). In the rate constant expressions, the barrier height is then ZPEexclusive. As stated above, this is the choice we make. (If, instead of using bottom of the potential well, one uses ZPE(s) of reactant(s) as the zero of energy, then the numerator of the vibrational partition function is 1, and the barrier height is ZPEinclusive, i.e., it is based on V^{G}_{a}.)
The quantized CVT rate constant (which we denote as CVT/T, where “/T” means tunneling) can be then written as

 (40) 

 (41) 
where
κ^{T} is the (temperaturedependent) tunneling transmission coefficient, which will be discussed in Section 2.3.
is the classical barrier height evaluated as the difference between the ground state energy of the variational transition state and the reactant(s), and because the latter is the zero of energy,
is simply the potential energy of the variational transition state;
Q^{GT} and
Q^{A} are the translationalmotionexcluded quantized partition functions;
Φ^{R} is the product of the reactants’ partition functions per volume, which is equal to

 (42) 
in which
m_{A}m_{B}/
m_{GT} is just the reduced mass of relative translational motion of A and B (because
m_{GT} =
m_{A} +
m_{B}). An alternative (equivalent) way for computing partition functions is that, in
eqn (40) and (41), the partition functions
Q^{GT} and
Q^{A} are the total partition function (including translational partition function), and
Φ^{R} is the product of reactants’ total partition functions per volume.
The above eqn (40) and (41) can also be rewritten as:

k^{CVT/T}(T) =κ^{T}Γ^{CVT}k^{TST}  (43) 
where
k^{TST} is the conventional TST rate constant, and
Γ^{CVT} is the canonical variational transition state theory recrossing transmission coefficient. (A pioneering physical picture for the transmission coefficients was developed in Hirschfelder and Wigner's work
^{245} in the very early days.) If we ignore tunneling (
i.e., we set
κ^{T} = 1), then the CVT variational transmission coefficient is

Γ^{CVT} = k^{CVT}/k^{TST}  (44) 
The value of Γ^{CVT} characterizes the importance of variational effects (recrossing) in the studied system, and it is smaller than unity. Under the harmonicoscillator approximation, CVT variational transmission coefficient is computed as

 (45) 
where
Q^{VTHO} is the vibrational partition function of variational transition state,
Q^{‡HO} is the vibrational partition function of the saddle point;
V^{VT} is the potential energy of the variational transition state, and
V^{‡} is the potential energy of the saddle point. The location of the GT is temperaturedependent, and so does the
Γ^{CVT}.
Similarly, for microcanonical variational transition state theory (μVT), we have

 (46) 
where
Γ^{μVT}(
E) is the microcanonical variational transmission coefficient. Microcanonical variational transition state theory provides a variational extension of the Rice–Ramsperger–Kassel–Marcus (RRKM) theory, which is based on conventional transition state theory. If quantum mechanical tunneling through rovibrational excited states of the transition state is considered, the number of states of the generalized transition state would be replaced by the cumulative reaction probability,
^{77,78} which thereby becomes an effective number of reactive states.
Notice that in the above equations, the overall rotational symmetry number of the molecule is included in the rotational partition function. Sometimes, the rotational symmetry numbers for TS and for reactant(s) in their rotational partition functions are factored out to the reaction symmetry number^{79}

 (47) 
and the corresponding VTST rate constants are expressed as (in which the rotational partition function excludes the rotational symmetry number)

 (48) 

 (49) 
where
σ^{R} and
σ^{‡} are the rotational symmetry number for the reactant and transition state structure, which are equal to the order of the rotational subgroup for polyatomic molecule; rotational symmetry number is 1 for
C_{∞v} point group (for instance, a heteronuclear diatomic molecule) and 2 for
D_{∞h} point group (for instance, a homonuclear diatomic molecule); for monoatomic species, which do not have rotational motion,
σ^{R} is taken to be 1.
A common misunderstanding about the reaction symmetry number is that it represents the number of equivalent ways (or chemically equivalent positions) for the reactants to react; this is simply not true, although in some cases they (the number of equivalent ways to react and the reaction symmetry number) might coincidentally have the same value. For instance, in the case of CH_{4} + OH → CH_{3} + H_{2}O, the number of equivalent ways for OH to abstract the H atom from CH_{4} is 4 (CH_{4} has four chemically equivalent H atoms), but the reaction symmetry number for forward reaction is 12 (the order of rotational subgroup for T_{d} point group (CH_{4}) is 12, and rotational symmetry number for OH is 1 and for TS (C_{s} point group) is 1); but in the case of CH_{4} + O → CH_{3} + OH, because the point group of the TS is C_{3v} (for which the order of rotational subgroup is 3), the forward reaction symmetry number is 4, which coincidently happens to be equal to the number of equivalent H atoms on CH_{4}. The reaction symmetry number should also include the contributions from the nonsuperimposable mirror images.^{80} However, if the MST method^{137,138} (which will be discussed in later section) has been applied, then the reaction symmetry number should exclude the mirror images, because they have already been included in MST partition function and we should not double count them.
We should remind ourselves that after this nonrigorous quantization of the classical canonical variational transition state theory, the CVT rate constant no longer necessarily provides an upper bound to the exact quantum rate constant. Nevertheless, the variationally optimized rate constants provide much better predictions than the ones obtained by putting dividing surface at the saddle point as in the conventional TST.
2.3 Quantum mechanical tunneling
2.3.1 Importance of tunneling.
Tunneling is the quantum mechanical process by which a particle has a finite probability to pass through a barrier, even though its energy is lower than the barrier. Nonclassical reflection refers to the quantum diffraction phenomenon that when the incident particle's energy is higher than the barrier, there is a finite probability of the particle being reflected. Because the treatment presented so far treats the reactioncoordinate motion as classical, it does not include either quantum mechanical tunneling through the reaction barrier or nonclassical reflection from the barrier. Tunneling increases the rate constant, and nonclassical reflection decreases it. Since tunneling occurs at lower energies where the Boltzmann factors are larger, tunneling is usually the more important one of the two phenomena for thermal rate constants, and the net effect of including these two quantum effects on reactioncoordinate motion are usually lumped together in a transmission coefficient that is simply called the tunneling transmission coefficient.
Quantum mechanical tunneling is almost always significant for reactions involving transfer of light atoms (such as hydrogen atoms, protons, or hydride ions) at low temperatures. It is often most easily noticed in large values of H/D kinetic isotope effects (KIEs). Sometimes it is also significant for heavieratom motion, for example, carbon tunneling.
2.3.2 Tunneling transmission coefficient.
The tunneling transmission coefficient κ(T) for VTST is defined as the ratio of the thermally averaged quantum tunneling probability to the quasiclassical transmission probability. The quasiclassical transmission probability is the one that is assumed in quasiclassical VTST, and it is a Heaviside step function that increases from 0 to 1 at the maximum of the effective potential. The tunneling transmission coefficient is: 
 (50) 
where E_{0} is the higher of either the groundvibrationalstate energy of the reactants or the groundvibrationalstate energy of the products (this is the lowest energy at which reaction can occur), P^{T}(E) is the quantum tunneling probability, and P^{C}(E) is the classical transmission probability. We assume that the effective potential for tunneling is the vibrationally adiabatic groundstate potential energy curve V^{G}_{a} of eqn (39), whose maximum value is called V^{AG}. Then, if the transition state were at the maximum of the vibrationally adiabatic groundstate potential energy curve, we would have 
 (51) 
where V^{AG} is the maximum of the vibrationally adiabatic groundstate potential energy V^{G}_{a} along the reaction coordinate. But in general, the CVT transition state is not at the maximum of the vibrationally adiabatic groundstate potential energy; thus we have to replace V^{AG} by where is the location of the CVT transition state at temperature T.
The above definition of the P^{C}(E) is consistent with the classical picture of a moving particle (i.e., no tunneling under the barrier and no reflection above the barrier). Carrying out the integral in the denominator of eqn (50), we get

 (52) 
As for the quantum tunneling probability P^{T}(E), it satisfies: 
 (53) 
The region below barrier is the tunneling region, and the region above barrier is nonclassical reflection region. The remaining problem is approximating the quantum tunneling probability P^{T}(E).
2.3.3 Semiclassical tunneling approximations.
In this section, we review various models for approximating the accurate quantum tunneling in reaction rate constant computations. Commonly used tunneling models can be divided to two classes: onedimensional and multidimensional tunneling models. The multidimensional tunneling models include the variation of the vibrational frequencies along the reaction path, as included by the final term of eqn (39). The onedimensional tunneling models do not. The onedimensional models reviewed in this article include the Eckart tunneling model,^{81,82} the Wigner tunneling model,^{83} the onedimensional truncated parabolic tunneling approximation;^{84} and multidimensional tunneling models include zerocurvature tunneling (ZCT),^{49,50} the parabolic model with an effective frequency,^{85} smallcurvature tunneling (SCT),^{86} largecurvature tunneling (LCT),^{71,87–90} microcanonically optimized multidimensional tunneling (μOMT),^{71} and leastaction tunneling (LAT).^{91–93} All of the onedimensional tunneling models and the ZCT tunneling model ignore the reactionpath curvature; this is not a good assumption for many cases. The effects of the reactionpath curvature often dramatically increase the amount of tunneling.
Tunneling models can also be divided into two classes in another way, namely models with analytic results (Wigner model, parabolic models, and Eckart model) and models where the tunneling probability must be calculated by numerical integration. The latter models are usually based on the semiclassical WKB approximation^{94} or – to be more precise – on a uniformized WKB approximation. The primitive WKB approximation suffers from infinities and has poor behavior for energies near the barrier top; uniformization^{95,96} is the process by which one include higherorder effects to obtain more uniformly accurate results, and when one does this, the results for tunneling through a barrier are often in good agreement (∼15% or better) with a quantal treatment of the same barrier model. But converting the multidimensional problem into a tractable problem with an effective barrier as a function of a tunneling coordinate (or a set of such tractable problems) is harder. We will begin with analytic onedimensional models.
2.3.3.1 Onedimensional tunneling models.
A popular approximation is to assume that the shape of the potential energy along the reaction coordinate s can be approximated by an Eckart potential:^{82} 
 (54) 
where 
y = exp[α(s − s_{0})]  (55) 

B = [V_{max}^{1/2} + (V_{max} − ΔV)^{1/2}]^{2}  (56) 
in which ΔV is the value of the potential at s = +∞ relative to s = −∞, V_{max} is the value of the potential at its maximum relative to s = −∞. And the constant s_{0} is defined to give the maximum value of V(s), (i.e., V_{max}) at s = 0: 
s_{0} = (2α)^{−1}ln(1 − ΔV/V_{max})  (57) 
The intrinsic barrier height V_{b} described by the potential is: 
V_{b} = V_{max} − max[0,ΔV]  (58) 
The parameters in Eckart model are usually chosen to fit the enthalpy of reaction at 0 K and the enthalpy of activation at 0 K. The parameter α determines the range of the potential, and it is related to the magnitude of the second derivative of the potential at its maximum by 
 (59) 
where the second derivative F_{s}, in practice, can only be approximately computed by: 
F_{s} = 4π^{2}μν^{‡}^{2}  (60) 
where ν^{‡} is the magnitude of the imaginary vibrational frequency of the saddle point, and μ is the reduced mass associated with the imaginaryfrequency vibrational mode. The reduced mass for reaction coordinate tunneling is often arbitrarily chosen to be 1 amu (atomic mass unit), or 1.008 amu for hydrogen atom tunneling (if we assume that the mass of the heavy atom is infinitely large). This way of approximating F_{s} is artificial, and it simply cannot be done unambiguously. The only advantages of the Eckart tunneling model are that one does not need to compute a reaction path, and the tunneling probability P^{T}(E) can be analytically evaluated, which provides computational simplicity, but it may yield quantitatively unsatisfactory results. The P^{T}(E) given by Eckart potential is: 
 (61) 
where 
a = 2π(2μE)^{1/2}/(ħα)  (62) 

b = 2π[2μ(E − ΔV)]^{1/2}/(ħα)  (63) 

d = 2π[2μB − ħ^{2}α^{2}/4]^{1/2}/(ħα)  (64) 
The Eckart tunneling model can be viewed as a crude approximation to the zerocurvature tunneling (ZCT) approximation, which systematically underestimates the tunneling due to neglecting reaction path curvature; however, the Eckart tunneling model often overestimates the ZCT tunneling transmission coefficient, which in turn partially compensates its intrinsic error. Eckart tunneling model uses a fit to V_{MEP}, in the sense that one fits it to enthalpy at 0 K for the energetics, but the frequency is still fit to the one based on V_{MEP}, which is probably why it often overestimates the tunneling probability compared to ZCT. (The V_{MEP} curve tends to be thinner for the nearly symmetric barriers that have the greatest tunneling.)
Although P^{T}(E) for Eckart tunneling can be expressed analytically, the tunneling transmission coefficient κ^{T}(T) needs to be computed by numerical integration. To further simplify the computation of κ^{T}(T), one could make additional assumptions (which may not be appropriate though) based on Eckart tunneling. If we consider a symmetric Eckart potential function, that is, we set ΔV = 0 in the asymmetric potential given by eqn (54) and we can get

 (65) 
The corresponding tunneling probability, which can be derived from eqn (61), is 
 (66) 
where 
 (68) 
If we assume the barrier is very high or very broad such that γ ≫ 1, then we have 
 (69) 
Plug the above tunneling probability into eqn (52), we can get 
 (70) 
in which 
x = γ(E − V_{max})/V_{max}  (71) 
Because of our previous assumption (γ ≫ 1), we carry out the integral in eqn (70) from −∞ to +∞, with a further assumption that the temperature T is very high,^{22}i.e., 
 (72) 
and finally an analytical expression can be arrived: 
 (73) 
This expression could at most be valid when the temperature T is very high and the barrier is very high or broad, but actually this simple expression for computing the tunneling transmission coefficient is basically never valid primarily because it uses a parabolic fit to V_{MEP} – not to V^{G}_{a}.
Eqn (73) can be further simplified since t/sin(t) ≈ 1 + t^{2}/6 for t ≪ 1, we obtain the simple Wigner tunneling approximation:^{83}

 (74) 
The Wigner tunneling transmission coefficient is very simple to compute; it only requires the imaginary frequency of the saddle point, but it is not useful for quantitative results.
It can be shown that, eqn (73) can also be derived from an infinite untruncated parabolic barrier.^{97} An untruncated parabola has E_{0} = −∞, and when T is lowered, the tunneling through an untruncated parabola become so unphysically large that the Boltzmann average in eqn (52) fails to converge and the transmission coefficient blows up when

hν^{‡}/(2k_{B}T) = π  (75) 
Skodje and Truhlar^{84} have eliminated this issue by using the uniformedized WKB semiclassical tunneling approximation, and they have further generalized the result to be applicable for any unsymmetrically truncated parabolictype barrier that is defined as: 
 (76) 
and the intrinsic barrier V_{b} is 
V_{b} = V_{max} − max[0,ΔV]  (77) 
The Skodje–Truhlar onedimensional truncated parabolic tunneling transmission coefficient κ^{T}(T) is given by: 
 (78) 
which is a singularityfree generalization of the result given by eqn (73). The above results (eqn (76)–(78)) are for onedimensional tunneling; later, Skodje et al. also showed that this parabolic model could also be generalized so that it includes the reactionpath curvature with the smallcurvature approximation,^{85} however, the correct effective mass for the SCT tunneling approximation was only derived later;^{86,98} it is not the one that is presented in this earlier paper^{85} (which is valid only if the reactionpath curvature is confined to a single transverse mode, which basically occurs only for atom–diatom reactions with collinear reaction paths). We will present the SCT tunneling model in the next section.
We do not recommend using any of the onedimensional models for quantitative studies; and they were reviewed in this subsection for completeness. Although these models are computationally inexpensive and certainly very easy to use, their accuracy is generally not good, and therefore for quantitative study, we recommend more accurate multidimensional tunneling models.
2.3.3.2 Multidimensional tunneling (MT) models.
In real dynamics, the motion of the reaction coordinate is not decoupled from the other degrees of freedom, and the curvature of the MEP in the massscaled coordinate system is one manifestation of the coupling.^{41,99} A major consequence of reactionpath curvature is that tunneling occurs in a corner cutting fashion.^{6,40,41,71,84–93,98–118} That is, the tunneling path does not need to follow the MEP;^{49} internal centrifugal effects cause the dynamical path to be displaced from the MEP; specifically, the quantum centrifugal effect makes the system move toward the concave side (that is the inside) of the MEP, not the convex side (that is the outside) as the classical centrifugal effect does. This is called the corner cutting effect, and it is due to the negative kinetic energy (imaginary momentum) in the classically forbidden region leading to a negative internal centrifugal effect.
First, we begin with the zerocurvature tunneling approximation (ZCT). In ZCT, the reaction path curvature is neglected, however, the contributions of other degrees of freedom are included in the vibrationally adiabatic potential V^{G}_{a}(s), which contains the vibrational groundstate energies of all the vibrational modes that are orthogonal to the reaction coordinate s, and this zeropoint energy depends on s. The contribution of vibrations transverse to the reaction path to V^{G}_{a}(s) is the reason why the SCT approximation is multidimensional. The ZCT tunneling probability is given by:

 (79) 
where E_{0} is the maximum between the V^{G}_{a} of the reactants and of the products; s_{*} is the location of the variational transition state, where the V^{G}_{a} reaches its maximum. The magnitude of the imaginary action integral is 
 (80) 
where s_{<} and s_{>} are the two turning points, at which V^{G}_{a} = E. The integrand is the imaginary part of the momentum of mass μ (that is the reduced mass of the isoinertial coordinate system). The ZCT approximation underestimates the tunneling transmission coefficient because it neglects reactionpath curvature; there are shorter tunneling paths than the MEP that have the same effective potential,^{119} which is given by eqn (39), and those paths dominate the tunneling.
In the smallcurvature tunneling approximation (SCT), an effective reduced mass μ^{eff} is introduced in computing the tunneling probability as expressed in eqn (79). The value of μ^{eff} is smaller than μ, which (as can be seen from eqn (80)) is equivalent to shortening the tunneling path. In the current version of SCT (which is also called the centrifugaldominant small curvature approximation), which is used in modern VTST calculations, μ^{eff} is determined by reactionpath curvature as follows:^{86,98}

 (81) 
in which ā(s) = κ(s)(s), and the reactionpath curvature κ(s) is , where B_{mF}(s) is the curvature component of vector B_{F}(s), which represents the coupling between the reaction coordinate s and a mode m that is perpendicular to it: 
 (82) 
where n_{i}(s) is the component i of the unit vector perpendicular to the generalized transition state dividing surface at s and L^{GT}_{i,m}(s) is the ith component of the eigenvector for vibrational mode m perpendicular to n(s) at s.
The effective harmonic potential V for the vibrational mode u_{1} (which is constructed by a local rotation of the vibrational axes such that B_{F}(s) lies along with it, and by construction, the curvature coupling in all other vibrational coordinates u_{2} to u_{F−1} are zero in this rotated coordinate system) is

 (83) 
where ω_{m}(s) is the vibrational frequency of mode m at reaction coordinate s. The associated turning point (s) for zeropoint motion in this harmonic potential is then: 
 (84) 
This form of turning point is able to avoid unphysically large tunneling probability that existed in the old version of the theory,^{85} which is caused by the assumption that all modes are extended to their turning points along the tunneling path; we no longer use the effective mass formula from the old version of the theory.
The turningpoint derivative terms dt_{m}(s)/ds are evaluated at s by a twopoint centraldifference method using the turningpoint results obtained at adjacent locations at which the normal modes are found. Because the gradient of the potential energy at the saddle point is zero, the effective reduced mass μ^{eff}(s = 0) is computed by a linear interpolation between the two closest values of μ^{eff}(s) as the following

 (85) 
where δs is the step size of the reaction coordinate.
The cornercutting model that led to the effective reduced mass μ^{eff}(s) does not apply to energies above the maximum of V^{G}_{a}, it is therefore inappropriate to compute the transmission probability at the nonclassical reflection region based on P^{SCT}(E). The transmission probability in the nonclassical reflection region is computed by

P^{SCT}(E) = 1 − P^{SCT}(2V^{G}_{a} − E), V^{G}_{a}(s_{*}) < E ≤ 2V^{G}_{a}(s_{*}) − E_{0}  (86) 
which means that the effective reduced mass is only used for computing the transmission probability at tunneling energies. SCT is currently the most widely used tunneling approximation among the many multidimensional tunneling models; and it is usually able to give a satisfactory tunneling transmission coefficient with a balanced computational cost (when compared to other more complex multidimensional models, which could be expensive for medium to large sized molecules).
Nevertheless, the tunneling path can be further optimized, and in some cases the system can tunnel directly into vibrationally excited states (i.e., it is not necessary for tunneling to be governed by the groundstate vibrationally adiabatic potential energy curve). The largecurvature tunneling approximation^{71,87–90} (LCT) was developed for more appropriately (as compared to SCT) describing reactions with large curvature of the reaction path. One of the key features in the LCT approximation is that the tunneling paths are set of straight lines that connect reactant and product sides, and these tunneling paths can connect vibrationalgroundstate turning points on the reactant side (in the exoergic direction of reaction) with vibrationally excited turning points on the product side. In the LCT method, the tunneling paths can be in a wider region than the nearly quadratic valley around the MEP, and the region that they pass through is called the reaction swath. The LCT tunneling probability includes the contribution from all the tunneling paths originating from the reactant valley and some of these occur earlier than when the system gets all the way in to where it hits the barrier. That is, in LCT, one can have a significant probability of tunneling before reaching the reactioncoordinate turning point. The classical reactioncoordinate turning points are the points at which the groundstate vibrationally adiabatic potential is equal to the total energy E in either the reactant and product valley. Unlike the treatment in ZCT and SCT, the LCT tunneling at energy E takes place all along the entrance channel up to the turning point, and it is initiated by the vibrational modes that are orthogonal to the reaction coordinate rather than the motion along the reaction coordinate. For large reactionpath curvature, the LCT approximation gives larger tunneling probability than SCT due to the very significant shortening of the tunneling paths; however, in the limit of small reactionpath curvature, the tunneling probability P^{LCT}(E) could be smaller than P^{SCT}(E), because SCT is a more accurate treatment in this limit region. At intermediate reactionpath curvature, SCT and LCT give similar results.
In the microcanonically optimized multidimensional tunneling model (μOMT), the tunneling probability at energy E is simply set to be the maximum between P^{LCT}(E) and P^{SCT}(E). When one uses μOMT, one often finds some contribution form largecurvature tunneling at low energy, even if the more important higherenergy tunneling is dominated by smallcurvature tunneling. In some older work, we used a method now called the canonically optimized multidimensional tunneling model (OMT), in which the thermally averaged tunneling transmission coefficient at temperature T is set to be the maximum between κ^{LCT}(T) and κ^{SCT}(T).
The most complete approximation for the quantum tunneling is the leastaction tunneling model (LAT).^{91–93} In this model, a sequence of paths that are parameterized by a parameter α is considered for each tunneling energy E and each final vibrational state (not necessarily the vibrational ground state). The parameter α is defined so that these LAT paths are all located at or between the MEP (α = 0) and LCT paths (α = 1). At every tunneling energy E, within the chosen sequence of paths, the tunneling path that has the largest tunneling probability is chosen to be the LAT tunneling path; that is, the LAT tunneling path is variationally optimized to minimize the action integral at E, with respect to the parameter α. The LAT tunneling was originally computationally very demanding, because one has to compute the information of huge number of points on potential energy surface (not just on MEP) including Hessians, but we have derived efficient algorithms for it.^{93} However, it is seldom used because the results are not sufficiently more accurate than those obtained by the μOMT method.
In the computations of the multidimensional tunneling transmission coefficients, it is very important for the computed reaction path (i.e., the range of s being computed) to be long enough so that the quantum tunneling is converged; otherwise a nonconverged quantum tunneling calculation can yield erroneous tunneling transmission coefficients. In practice, one can extend the range of the reaction path that was computed in the previous calculation and then redo the computation (with some possible restart options to take advantage of the information produced in the previous run), and the tunneling transmission coefficients should remain the same at all temperatures. Extrapolation techniques are available to minimize the needed length of reaction path computation.^{120}
2.3.4 Tunneling and the treatment of reaction intermediates.
An important issue involving tunneling is the treatment of reaction intermediates. Consider a twostep reaction: 
 (87) 
The first step is the formation of an intermediate (I), with a forward rate constant of k_{1}; the intermediate I further proceeds to yield the product P, with a rate constant of k_{2}. The overall reaction rate for the formation of product is: 
 (88) 
and if the intermediate I is in equilibrium with A and B, we would have 
 (89) 
where K is the equilibrium constant for the first step, and K = k_{1}/k_{−1}. And thus the overall bimolecular rate constant isWith (variational) transition state theory, the rate constant k_{2} is 
 (91) 
where κ^{T}_{2} is the tunneling transmission coefficient computed by using the intermediate I as the reactant. The equilibrium constant K of the first step is 
 (92) 
and thus, the overall bimolecular rate constant is equivalent to 
 (93) 
If one ignores the intermediate, that is if one computes the rate constant directly for A + B → P, the rate constant expression is almost the same as eqn (93); the difference is that κ^{T}_{2} in eqn (93) is replaced by κ^{T}_{1}, which is the tunneling transmission coefficient computed by using A and B as the reactants. To discuss this, we consider the case of an exothermic reaction at low temperature where the reactant is in its ground state. If the top of the effective barrier for tunneling is below the zero point energy of the reactants, then there is no tunneling because all systems have enough energy to pass over the barrier without tunneling. Thus κ^{T}_{1} is unity (or less than unity if we include nonclassical reflection). But if the pressure is large enough to stabilize states of the intermediate that have energies lower than the reactant zero point energy and lower than the top of the effective barrier for tunneling, those low energy states can tunnel and κ^{T}_{2} can be greater than unity. This shows that the wellknown dictum that one can ignore an intermediate that occurs before the ratedetermining transition state is only partially true, and it is dangerous to ignore the existence of a van der Waals intermediate in a bimolecular reaction (or similarly, in a twostep unimolecular reaction). This is a very general issue since most bimolecular reactions have barrier, and a groundelectronicstate barrier is always preceded by a van der Waals minimum. However, the issue is usually important only for reactions with low barriers. For a high barrier, since the tunneling usually has important contributions only from energies at most a few kcal mol^{−1} below the barrier top, the tunneling from intermediate states whose energy is below the reactant zero point energy is expected to be negligible.
There is another complication, namely the assumption that the intermediate is in equilibrium with A and B. If one is not in the highpressure limit, there may be inadequate collisional stabilization of the complex. In the lowpressure limit, there are no stabilizing collisions, and therefore there are no tunneling contributions at energies below the groundstate energy of the reactants^{121} (which would be the lowest possible energy that a complex could have). The lowpressure limit and the highpressure limit are easily handled as limiting cases, but at intermediate pressures one must use the theory of pressuredependent reactions, which is treated in a later section.
2.4 Recrossing transmission coefficient
The basic assumption of quasiclassical transition state theory, which may be used to derive the equations in Section 2.2, is that the systems with energy higher than the vibrationally adiabatic potential barrier (which may be called the quasiclassical threshold energy) react with unit probability and systems with lower energy do not react, i.e., the reaction probability is a Heaviside step function as a function of energy in the reaction coordinate. Because the vibrationally adiabatic barrier involves vibrational frequencies and rotational moments of inertia that depend on the reaction coordinate, VTST goes beyond the assumption of a separable reaction coordinate that is assumed in purely classical conventional transition state theory.
In Section 2.3, we saw the use of a tunneling transmission coefficient to account the breakdown of the step function probability due motion along the reaction coordinate being quantum mechanical so there is a small probability that systems with higher energy than the effective barrier diffract off the barrier top and do not react and a small probability that systems with less energy that the barrier tunnel through it. The latter events happen at lower energy where the Boltzmann factor is higher, and so their effect usually dominates over the effect of nonclassical reflection, and the tunneling transmission coefficient is usually great than unity. We take account of tunneling by doing a multidimensional semiclassical calculation that accounts for two aspects of reactioncoordinate nonseparability, namely (i) the alreadymentioned reactioncoordinate dependence of the vibration frequencies and (ii) the corner cutting effect by which the tunneling path involves displacement of the transverse vibrations to the concave side of the curved minimum energy path during the tunneling event.
There is another reason why the step function probability assumption can break down, namely that even for classical reaction coordinate motion, the nonseparability of the reaction path can lead to less than unit probability of reaction for systems with energy greater than the quasiclassical threshold because systems cross the transition state but then recross it, making the net reactive flux smaller than the oneway flux. We can account for this by including another transmission coefficient in the rate constant; we call this the recrossing transmission coefficient. Note that the recrossing transmission coefficient could also be called the quasiclassical transmission coefficient because the effect is present even for classical reactioncoordinate motion.
In general, then there are two contributions to the transmission coefficient, and we write
where
γ is the full transmission coefficient,
κ is the tunneling transmission coefficient, and
Γ is the recrossing transmission coefficient. We can apply this to the quasiclassical VTST rate constant calculated without tunneling or recrossing to get a better answer. When we do this, we usually set
Γ = 1. For example, using canonical variational transition state theory, we write the rate constant with tunneling as
where
k^{CVT} is the quasiclassical VTST rate constant. However, sometimes we use a different notation, namely

k^{CVT/T} = κΓ^{CVT}k^{TST} = γk^{TST}  (96) 
where
k^{TST} is the quasiclassical conventional VTST rate constant, and

Γ^{CVT} = k^{CVT}/k^{TST}  (97) 
These different conventions should be clear from the context.
The quantity Γ^{CVT} is sometimes called the variational effect because it accounts for the difference between variational transition state theory and conventional transition state theory when both are applied without tunneling. Its deviation from unity is a measure of the extent to which VTST accounts for the breakdown of the norecrossing assumption of quasiclassical transition state theory at the conventional transition state that passes through the saddle point. For example, a recrossing transmission coefficient of 0.8 means that 20% of the quasiclassical flux through the conventional transition state does not contribute to the net oneway flux because the true dynamical bottleneck is the variational transition state, not the conventional transition state. This would mean that 20% of the flux through the conventional transition state is nonreactive because it does not pass through the variational transition state, which has a higher free energy of activation.
One can also define an exact recrossing transmission coefficient as

Γ^{exact} or γ^{exact} = k^{exact}/k^{TST}  (98) 
This is the definition implicitly used in original literature,
^{3} where it is simply called the transmission coefficient.
2.5 Importance of using curvilinear coordinates in VTST/MT calculations for vibrations transverse to the reaction coordinate
In the VTST/MT calculations, the computations of both the variational transition states and the multidimensional quantum tunneling transmission coefficients rely on the computed vibrationally adiabatic groundstate energy along the reaction path; the generalized normal mode analysis is performed at evenly spaced points along the minimum energy path. For nonstationary points (notice that the reactants, products, and the conventional transition state structures are stationary points) along the reaction path, the frequency calculations have to be carried out based on projected Hessians (which project out the translational modes, overall rotational modes and the mode corresponding to the reactioncoordinate degree of freedom). One should keep in mind the fact that the projected Hessians in this generalized normal mode analysis are coordinatesystemdependent,^{75,76} and so are the boundmode vibrational frequencies of nonstationary points.
The reason that the vibrational frequencies along the MEP are coordinatedependent can be explained as follows. The choice between the rectilinear coordinates (which are linear combinations of Cartesian coordinates) and the curvilinear coordinates (internal coordinates) is equivalent to the choice between two different definitions, s and s′, of the reaction coordinate for points off the MEP; they have the same values for stationary points on the MEP, but not for nonstationary points. These two different choices of the reaction coordinate are related by^{76,122}

 (99) 
in which
F is the number of internal coordinates, and
F − 1 excludes the reaction coordinate
s (which is also called
q_{F}). The
q_{i} represent curvilinear (internal) coordinates, which measure distortions away from the MEP, and by definition they vanish on the MEP. The
b_{ij} involve secondorder partial derivatives of
s′ with respect to
q_{i}, and they are in general nonzero. The double summation runs over a complete, nonredundant set
q_{1},
q_{2},…,
q_{F−1} of such coordinates that are orthogonal to the MEP. The matrix elements of the Hessians expressed using these two different reaction coordinate definitions are related by:
^{123} 
 (100) 
in which
q ≡ {
q_{1},
q_{2},…,
q_{F−1},
s} and
q′ ≡ {
q_{1},
q_{2},…,
q_{F−1},
s′}. Clearly, at stationary points, where the firstorder derivative of
V with respect to
s is zero, the two Hessians are the same; Cartesian coordinates and internal coordinates yield the same frequencies for the stationary points. However, at nonstationary points, where the (∂
V/∂
s) term does not vanish, the frequencies are coordinatesystemdependent.
And thus, the important consequence is that, the results of the VTST/MT calculations depend on the choice of the coordinates^{75,76} that are used for the generalized mode analysis along MEP. It is highly possible for a set of Cartesian coordinates (a linear combination of which yields a set of rectilinear coordinates) to unphysically yield one or more imaginary frequencies along the reaction path^{75,76} (along the reaction path, in the absence of a real bifurcation, there should be no imaginary frequency after the projection); and the resulting VTST/MT rate constants can be significantly affected, and thus they could be unreliable or even wrong. Along the reaction path, the free energy cannot be computed at the regions where one has imaginary frequencies, and thus Polyrate program^{42} arbitrarily sets the vibrational partition function to very large number (which decreases the Gibbs free energy by a large amount) in order to ensure that the computed canonical variational transition state will never occur in these unphysical regions which contains imaginary frequencies. Using Cartesian coordinates for generalized normal mode analysis along the reaction path could possibly cause: (1) unphysical discontinuity of the vibrationally adiabatic groundstate energy (V^{G}_{a}) along the reaction path; (2) unphysical sudden jump(s) on the Gibbs free energy at a given temperature along the reaction path, because of the abovementioned arbitrary setting that the program has to make; (3) in the case of microcanonical VTST calculations, using Cartesian coordinate can cause unphysical sudden jump(s) in the computed vibrational number of states N^{μVT}(E,s) and thus the energydependent microcanonical variational transition state would be located in these unphysical regions, which leads to unphysical results for the μVT rate constants.
In order to obtain reliable VTST/MT rate constants, one should first choose a set of physically reasonable curvilinear coordinates that are capable of describing the atomic motion along the whole reaction path, and then use this set of coordinates to carry out the generalized normal mode analysis and to further compute the variational transition states and multidimensional tunneling. In Polyrate, one could choose either a set of 3N − 6 (where N is total number of atoms) nonredundant internal coordinates^{75} (via the “CURV2” option), or redundant internal coordinates^{124} (via the “CURV3” option), as curvilinear coordinates for performing the generalized normal mode analysis. A set of physically meaningful and wellchosen internal coordinates is usually able to give an imaginaryfrequencyfree reaction path in physically important regions (i.e., the regions where the variational transition states are located and where the quantum tunneling take place). In some cases, it is useful to also reorient the generalizedtransitionstate dividing surface (via the “RODS” option) to eliminate the imaginary frequencies along the path.^{125} As for the computations of the multidimensional tunneling transmission coefficients, it is very important for the computed reaction path (i.e., the range of s being computed) to be long enough so that the quantum tunneling is converged; otherwise a nonconverged quantum tunneling calculation can yield erroneously large tunneling transmission coefficients. In practice, one can extend the range of the reaction path that was computed in the previous calculation and then redo the computation (with some possible restart options to take advantage of the information produced in the previous run), and the tunneling transmission coefficients should remain the same at all temperatures.
2.6 Stateselected rate constants
Although transition state theory is formulated for the calculation of the reaction rates of thermally equilibrated or microcanonically equilibrated reagents, extensions have been proposed to approximate vibrationalstateselected rate constants.^{126,127} Liu recently reviewed the experimental progress made in studying the modeselected reaction dynamics.^{128} In this section, we briefly review using VTST to study stateselected reaction rate constants.
There are two basic ways to extend VTST to reactions in which the initial vibrational state is specified to have some specific excitation. One is to assume vibrational adiabaticity of the excited mode;^{127} then, in the sums over vibrational quantum numbers for the stateselected mode, only one term contributes. Another is based on nonadiabaticity. We first review the former approach and then the latter one.
Consider the vibrational partition function with F vibrational degrees of freedom assumed separable (for stationary points, F includes all the vibrational degrees of freedom; and for nonstationary points, F includes all the boundmode vibrational degrees of freedom that are orthogonal to MEP, which excludes the reaction coordinate),

 (101) 
where
m labels a vibrational mode, and
n_{m} is the number of quanta in that mode, and
n^{max}_{m} is the number of quanta used to converge the summation. If the vibrational energy levels are computed by harmonic oscillator (HO) approximation

 (102) 
where
v_{m} is the vibrational frequency of mode
m. This yields the HO vibrational partition function of
eqn (37), which is used in computing thermal VTST rate constants in the harmonic approximation. In modeselected VTST, one modifies
eqn (101) to

 (103) 
in which the modes from
m = 1 to
m =
F^{r} are the selected modes (
i.e., each of these modes is restricted to have a specific number of excitation quanta equal to
n^{r}_{m}); and the remaining
F–
F^{r} modes are treated as usual with the summation over all the vibrational quantum numbers
n_{m′} of modes
m′ going up to a high enough value to converge the partition function,
i.e., up to a large number that converges the summation. This restricted vibrational partition function is a function of the selected vibrational quantum states, which are characterized by the set of the specified vibrational quantum numbers {
n^{r}_{m}}. In adiabatic modeselected VTST, the locations of the variational transition states are determined based on the vibrational partition functions of the generalized transition state using
eqn (103) with the boundmode vibrational frequencies along the reaction path. This same restriction is also used to generate stateselected conventional transition state theory, in which the generalized transition state is chosen to be at the saddle point rather than variationally optimized.
The multidimensional quantum tunneling transmission coefficients in stateselected VTST is computed based on the following stateselected vibrationally adiabatic potential:

 (104) 
The effective potential of eqn (104) is based on the adiabatic approximation, which is a good approximation when the dynamical bottleneck occurs prior to the region of large reactionpath curvature; large reactionpath curvature causes strong vibrational nonadiabaticity,^{129} and in such a case, one should use the second way to treat stateselected reactions, namely the diabatic one^{126} with the vibrationally diabatic potential given by

 (105) 
which is similar in form to the expression for the vibrationally adiabatic potential
V^{g}_{a}, but instead of using adiabatic frequencies, it uses diabatic frequencies which change smoothly along the reaction coordinate. When the adiabatic frequencies
v_{m}(
s) show an avoided crossing in the frequencies
v_{m}vs. s diagram (this is called the vibrationalmode correlation diagram, in which the curves corresponding to modes of the same symmetry do not cross), it is often a better approximation physically to assume that the modes preserve their character in passing through the avoided crossing region. For instance, in the OH + H
_{2} → H
_{2}O + H reaction,
^{126} there are five bound vibrational modes along the reaction coordinate
s, in which the three modes with the lowest adiabatic vibrational frequencies (of which the modes are denoted as
m = 1, 2, 3) change smoothly along
s, and the two highestvibrationalfrequency modes (which are denoted as
m = 4, 5) with the same symmetry have an avoided crossing at
s_{0} (which is about −0.16 Å). The diabatic frequencies
v^{d}_{m}(
s) are then defined as:

 (106) 
where
θ is the step function, and
δ is Kronecker delta.
Stateselected VTST calculations have been tested for various systems,^{130–136} and good agreement with experimental values or accurate quantum dynamics are confirmed. The most sophisticated treatment^{131} assumes that the reaction is vibrationally adiabatic up to the first occurrence in proceeding from reactants to products of an appreciable local maximum in the reactionpath curvature; at that point the reaction is treated as if all flux is suddenly diverted to the ground vibrational channel.
3. Multistructural and multipath variational transition state theory
3.1 Multistructural anharmonicity
The textbook way of computing molecular partition function is based on a single lowestenergy structure and on the (quasi)harmonic oscillator approximation. A harmonic oscillator has only one minimumenergy structure, but real species often have multiple structures, e.g., gauche and trans 1,2dichloroethane. This section discusses how to include those features in VTST. Some of the material in this section applies to any reactant or transition state that has multiple structures, but most of the discussion is written in the framework of torsional multistructural anharmonicity, which has two components: multiplestructure anharmonicity and torsional potential anharmonicity. The multiplestructure and torsional potential anharmonicity effects can be very important in the computation of partition functions, especially at high temperatures.
If a reagent or a transition state being studied has one or more torsional degrees of freedom (internal rotations), the harmonic oscillator approximation fails because (a) the potential function is not quadratic and (b) torsions usually generate multiple structures, each of which contributes to the partition function. Issue (a) is called torsional potential anharmonicity,^{137} and issue (b) is called multiple structure anharmonicity. Multistructural torsion anharmonicity can be treated with the recently developed multistructural torsion (MST) method, which has two versions: MST with uncoupled torsional potential method (MST(U) method), and MST with coupled torsional potential method (MST(C) method). We recommend the second version of this method, which is based on a coupled torsional potential.^{138} Including the MST treatment of anharmonicity in VTST will be called multistructural VTST (MSVTST) or multipath VTST (MPVTST), depending on how completely it is included.
A methyl group can usually be safely treated as nearly separable from other torsions, but more generally the torsional degrees of freedom in a molecule are coupled; they interact with each other, and sometimes this coupling is strong. If we rotate each torsional bond to n different angular positions, and if we have t torsions (temporarily excluding methyl groups because rotating methyl groups does not generate new distinguishable structures; however we do consider all the torsions including methyl groups in the final MST partition function), we will generate n^{t} initial structures. However, after optimizing all these initial structures by some electronic structure method (this step is called the conformational search), the distinguishable conformers one finds will not be equal to n^{t}. All of the distinguishable conformers that are generated by torsions are included in the computation of the MST partition function, including nonsuperimposable mirror images. The multiple structure effect can be very large when one is treating a long chain molecule.
One may propose to treat all the local minima (or in the case of transition state, all the first order saddle points on reaction paths from reactants to products) with the harmonic oscillator formulas; and this is called the multistructural local harmonic oscillator approximation (MSLH).^{139} However, treating lowfrequency modes that are torsions or have large contributions from torsions by the harmonic oscillator approximation can be very inaccurate. The harmonic potential is a good approximation to the real potential only around the equilibrium structure. The harmonic vibrational partition function has the hightemperature limit:

 (107) 
and this is the vibrational partition function of a classical harmonic oscillator.
At very low temperatures where k_{B}T ≪ hν, an internal rotation can be reasonably well treated with the harmonicoscillator approximation at each minimum of the torsional potential; at very high temperature where k_{B}T is much larger than the internal rotational barrier W, and the internal rotation becomes a free rotor. At intermediate temperature where hν < k_{B}T < W, torsional degrees of freedom should be treated as hindered rotors.
Here we review the MST method based on a coupled torsional potential (i.e., MST(C) method). By “coupled”, we mean that we are using consistently coupled torsional vibrational frequencies, coupled effective torsional barriers, and coupled internal moments of inertia; the kinetic couplings between the torsional motions and the overall rotation of the molecule are also fully treated.
The MST rotational–vibrational (which can be abbreviated as rovibrational) partition function is computed as the summation of the singlestructure torsional rovibrational partition function (SST) of the unique structure j. In the equations used for the MST partition function, a unique structure refers to a quantum mechanically distinguishable (both structurally and energetically) conformational structure, which is generated by a change in vibrational coordinates, usually by a torsion; the unique structures are torsional conformers – not chemical isomers. And a distinguishable nonsuperimposable mirror image that can be generated via torsion also counts as a unique structure. The MST partition function is:

 (108) 
where the SST rovibrational partition function for distinguishable structure
j is defined as

 (109) 
where
Q^{rot}_{j} is the rotational partition function of structure
j,
Q^{HO}_{j} is the harmonicoscillator (or quasiharmonicoscillator) vibrational partition function of structure
j,
U_{j} is the relative Born–Oppenheimer equilibrium energy of structure
j with respect to the Born–Oppenheimer equilibrium energy of the lowestenergy conformer (which is denoted as the global minimum or GM),
t is the total number of torsions, and
f_{j,η} is the torsional correction factor for the coupled torsion
η of structure
j. If we set all the torsional correction factors in
eqn (109) to unity, then this equation reduces to the singlestructure harmonicoscillator rovibrational partition function (SSHO) of structure
j.
The rotational partition function for structure j is computed using the classical rotor approximation, which is sufficient for most purposes, and which yields:

 (110) 
where
σ_{rot,j} is the overall rotational symmetry number for structure
j (each conformer
j has its own symmetry number), and
I^{rot}_{j,1},
I^{rot}_{j,2}, and
I^{rot}_{j,3} are the three principal moment of inertia associated with the overall rotation of the rigid structure
j. The rotational symmetry number of conformer
j is equal to the order of the rotational subgroup for a polyatomic molecule (conformer
j); the rotational symmetry number is 1 for a heteronuclear diatomic molecule and 2 for a homonuclear diatomic molecule.
The torsional correction factors of all the coupled torsions are treated together by

 (111) 
where
q^{RC(C)}_{j,η} is the reference classical partition function of coupled torsional motion
η for structure
j using the coupled effective torsional barrier, and
q^{CHO(C)}_{j,η} is the classical harmonic oscillator partition function using consistently coupled frequencies.
The reference classical partition function of structure j is an approximation to the fully coupled classical partition function Q^{FCC} (where all t torsions and the overall rotation are coupled):

 (112) 
where
σ_{rot} is the overall rotational symmetry number,
σ_{η} is the torsional symmetry number for torsion
η (the symmetry number for a torsion can be determined by treating the least symmetric of the two rotating fragments as a fixed frame and counting the number of identical structures obtained when the more symmetrical top is rotated from 0 to 360 degrees with respect to this torsional bond and relaxed),
V is the potential energy as function of torsional angles
ϕ_{η}, and det(
S) is the determinant of the Kilpatrick and Pitzer moment of inertia matrix for overall and internal rotation,
^{140} which is uniquely defined by the torsional coordinates (dihedral angles). The fully coupled classical partition function cannot be analytically integrated, and the computation of the whole potential energy
V as function of torsional angles is generally unaffordable for systems with multiple torsions, and therefore several approximations are made in order to evaluate
Q^{FCC}: (a) we approximate the determinant of
S(
ϕ_{1},…,
ϕ_{t}) matrix in the vicinity of each conformational structure
j by its value at a local minimum or the saddle point; (b) we use a congruent transformation to block diagonalize the
S matrix so that it can be written as the product of the three principal moment of inertia for the overall rotation of a rigid structure, with the torsional moment of inertial matrix
D (which fully accounts for kinetic coupling of torsions to each other and to overall rotation); (c) we approximate the potential
V(
ϕ_{1},…,
ϕ_{t}) with the uncoupled local reference potential in the vicinity of structure
j as

 (113) 
where, in the coupled MST calculation,
W_{j,τ} is the coupled effective torsional barrier
W^{(C)}_{j,τ} for torsion
τ, which is computed by diagonalizing the torsional Wilson–Decius–Cross
GF submatrix,
^{46} and
M_{j,τ} is the local periodicity of torsion
τ for conformer
j, which is related to the number of structures (including both distinguishable and indistinguishable ones) that can be generated by torsion
τ.
The local periodicity of an uncoupled torsional degree of freedom is the same in all unique conformers, and it is an integer equal to the number of distinguishable conformers (including mirror images) that can be generated by this torsion, multiplied by the torsional symmetry number of this torsion. For instance, for a nearlyseparable uncoupled –CH_{3} group, there is only one distinguishable conformer that can be generated by rotating this torsional bond, and the torsional symmetry number of methyl group is 3; and thus it has M = 3. In strongly coupled systems, M_{j,τ} is generally a conformerdependent noninteger, and it is assigned by using Voronoi tessellation.^{141–143} As a result of putting all these elements together, our reference classical partition function for structure j, in which all the coupled torsions are treated together in product form, is

 (114) 
where
I_{0} is modified Bessel function. Notice that, for coupled torsions, it is the product of local periodicities that appears here, rather than a separate dependence on the local periodicity of each coupled torsion.
When torsions are strongly coupled together, it is sometimes impossible to assign integer M_{j,τ} values for each torsion. We divide the t torsions into two types: nearly separable (NS, for instance CH_{3} groups) and strongly coupled (SC). We use the notation NS:SC = t_{NS}:t_{SC} to denote that t_{NS} torsions are treated as nearly separable and t_{SC} torsions are treated as strongly coupled. In general, the strongly coupled coordinates may be further partitioned into two or more subspaces, with each subspace involving only those coordinates that are strongly coupled to each other. Each of the strongly coupled subspaces is treated by Voronoi tessellation separately. Voronoi tessellation divides a space into cells around a discrete set of points, and the space to be tessellated here is described by the dihedral angles ϕ_{1}, ϕ_{2},…, ϕ_{tSC}, and the points correspond to structures. Each cell corresponds to a specific structure and consists of all torsional configurations closer to this structure than to any other structure when only the t_{SC} strongly coupled degrees of freedom are considered. In Voronoi tessellation, all the distinguishable and indistinguishable structures that can be interconverted via torsions are included. The volume of the SC torsional subspace is then given by

 (115) 
where
t_{SC} is the number of torsions in the strongly coupled space, and
σ_{τ} is the torsional symmetry number for torsion
τ. Whenever we use Voronoi tessellation, we assume that the torsional subspace is so strongly coupled that we cannot assign
M_{j,τSC} by considering each torsion separately; then we replace all
M_{j,τSC} for strongly coupled torsions of a given conformer
j by a single
M^{SC}_{j} equal to

 (116) 
As we can see in
eqn (111) and (114), it is the product of the local periodicities rather than the individual terms that are needed in the computation of the torsional correction. In the practical Voronoi algorithm,
^{141}eqn (116) is computed by:

 (117) 
In this sampling algorithm, a total number of
N_{tot} random points in the torsional space are generated, and each point is assigned to the structure that is the closest to it (in other words, each point is assigned to a particular Voronoi tile), and
N_{j} is the number of points that are assigned to structure
j.
Here we give some numerical examples of the M^{SC}_{j} values for the strongly coupled (SC) torsional degrees of freedom in some systems; all the methyl groups in the following examples are excluded in the given M^{SC}_{j} values, because they are treated as nearly separable with M = 3. For 1pentyl radical (for which the conformational search is carried out by M062X/MG3S), the M^{SC}_{j} values for the SC torsions in its eight distinguishable conformers (which excludes 7 distinguishable mirror images) are, from lowestenergy conformer to highestenergy conformer, respectively 2.67, 2.63, 2.73, 3.84, 3.82, 3.82, 3.36, and 3.39. The transition state structure of the hydrogen abstraction reaction from a methyl group of the tertbutanol by HO_{2} radical (the products are (CH_{3})_{2}C(OH)CH_{2}˙ radical and H_{2}O_{2}) possesses five SC torsions and has 46 distinguishable conformers (i.e., 23 pairs of nonsuperimposable mirror images, which are found by M08HX/MG3S); the M^{SC}_{j} values range from 2.00 to 2.45, with a mean value of 2.17 and a standard deviation of 0.13. For (S)secbutanol (for which the conformational search is carried out by M08HX/MG3S), the M^{SC}_{j} values for the SC torsions in its nine distinguishable conformers (from the lowestenergy conformer to the highestenergy one) are respectively 3.01, 2.97, 3.01, 3.00, 3.00, 2.95, 3.02, 3.05, and 3.00.
The classical harmonic oscillator partition function q^{CHO(C)}_{j,η} with coupled frequencies is:

 (118) 
where
F is the number of vibrational degrees of freedom,
ω_{j,m} is the vibrational frequency of normal mode
m, and
_{j,m} are the nontorsional vibrational frequencies; these frequencies are computed by diagonalizing respectively the full Wilson–Decius–Cross
GF matrix and its nontorsional submatrix with nonredundant internal coordinates.
The final result is

 (119) 
In general, for a system with multiple coupled torsional degrees of freedom, the accurate partition function is unaffordable to compute. Nevertheless, for twodimensional systems, i.e., molecules with two coupled torsional degrees of freedom, where the highly accurate partition functions are available via the extended twodimensional torsion method^{144} (E2DT) which is a method that is based on the solution of the torsional Schrödinger equation and that includes full coupling in both the kinetic and potential energy, the satisfactory accuracy of the MST method with coupled torsional potential, i.e., MST(C) method, has been confirmed; and it is concluded that “If the E2DT method is not affordable, the MST(C) method is the best option to calculate rovibrational partition functions including torsional anharmonicity”.
To carry out MST calculations, one has to do an exhaustive conformational search at the first step, which might be unaffordable for structures with a large number of torsional degrees of freedom, even when using economical density functional methods. Recently we have proposed a much less expensive method called duallevel MST,^{145} which requires the user to carry out a conformational search with a chosen lowlevel theory (for instance, a semiempirical molecular orbital method) and then reoptimize the 30% energetically lowest lowlevel structures at a higher level. Duallevel MST is able to combine the information from a lowlevel theory with that from a higherlevel theory to obtain an approximate MST partition function that agrees with the full highlevel MST partition function within a factor 2 to 3; an even more encouraging feature is that we have shown that when one is using duallevel MST partition functions in computing rate constants or thermodynamic functions (in which cases one usually needs the ratio of the partition functions rather than an individual partition function, and the errors of the duallevel method canceled out to a large extent), very satisfactory results (as compared to using full highlevel MST partition functions) can be obtained.
3.2 Multistructural variational transition state theory
In multistructural variational transition state theory^{146,150} (MSVTST), we replace the single structure harmonic oscillator partition functions in the expression of the VTST rate constant, by MST partition functions. And therefore, the multistructural and torsional anharmonicity is incorporated in the VTST theory.
We now derive the MSVTST rate constant for unimolecular reactions, for bimolecular reactions one can follow the similar approach. For unimolecular reactions, the single structural variational transition state theory rate constant with multidimensional tunneling can be written as

 (120) 
where
Q^{GTSSHO} and
Q^{RSSHO} are respectively the singlestructural (lowestenergy structure) harmonicoscillator rovibrational partition function for the generalized transition state and for the reactant. The location of the variational transition state is denoted as
s_{*};
V(
s =
s_{*}) is the potential energy of the variational transition state. We can rewrite this expression in terms of the recrossing transmission coefficient and the conventional (singlestructural) TST rate constant

 (121) 
where
Γ^{VTST} is the VTST recrossing transmission coefficient (which could be computed from canonical, microcanonical, or any other version of variational transition state theory) resulting from variational effects. In CVT, the recrossing transmission coefficient is a special case of
Γ^{VTST} and is computed as

 (122) 
where
Q^{‡SSHO} is the singlestructural harmonicoscillator rovibrational partition function for the lowestenergy saddle point. In order to take the MST effects into VTST calculation, we replace the singlestructural transition state theory rate constant with the multistructural TST (MSTST) rate constant in
eqn (121). The multistructural TST rate constant is defined by replacing the singlestructural partition functions with MST partition functions

 (123) 
which can be written in terms of the multistructural and torsional anharmonicity factor
F^{MST}_{fwd} of forward reaction:

k^{MSTST}(T) = F^{MST}_{fwd}k^{SSTST}(T)  (124) 

 (125) 
The multistructural canonical variational transition state theory (MSCVT) rate constant is defined as the following:

k^{MSCVT/T}(T) = κ^{T}Γ^{CVT}k^{MSTST}(T) = F^{MST}_{fwd}k^{SSCVT/T}(T)  (126) 
where the singlestructural CVT rate constant is

k^{SSCVT/T}(T) = κ^{T}Γ^{CVT}k^{SSTST}(T)  (127) 
Notice that, in the MSCVT theory presented here, we approximate the ratio of the MST partition function to the SSHO partition function at the variational transition state as the ratio at the conventional transition state; and the variational transmission coefficients are computed based on singlestructural (quasi) harmonicoscillator approximation. This is usually a very good approximation because usually the variational transition states are very close to the saddle point (at generally accessible temperatures, the CVT variational transition state is usually not more than 0.1 to 0.2 Å away from the saddle point). Also, by adopting this approximation, we completely eliminate the need for computing the multiplestructural anharmonicity (which requires optimizations and generalized normal mode analysis for all the possible generalized conformers of the nonstationary points) along the reaction path, which is simply unaffordable for large or even mediumsized transition state structures. For bimolecular reactions, the factor
F^{MST}_{fwd} for forward reaction is defined as:

 (128) 
For the reverse reaction, we can also define the F^{MST}_{rev} factor for the reverse reaction, in which one replaces the reactant partition function in the denominator with the product partition function. We can also define the MST standardstate reaction equilibrium constant K°^{,MST} being the product of F^{MST}_{rxn} factor for the reaction (which is equal to F^{MST}_{fwd}/F^{MST}_{rev}), with the singlestructural harmonicoscillator standardstate reaction equilibrium constant K°^{,SSHO} (which can be derived from SSHO standard Gibbs free energy of reaction ); and the MST standard Gibbs free energy of reaction .
The basic assumption of transition state theory is that the conformational structures of the reactants in phase space are in local equilibrium with each other and with those of the transition state that correspond to motion toward the product. MSVTST is applied under this assumption, and it assumes that the interconversion between the conformers of the reactant is much faster than the chemical reaction itself.^{121} However, if the interconversion is comparable or even slower than the reaction being considered, as possibly in the case where the torsional barriers are high as compared to the reaction barrier, then these different conformers should be considered as different reactants. It is also interesting to point out that, as compared to the generalized Winstein–Holness (GWH) equation,^{147,148} or other similar treatments (such as the socalled multiconformer TST), which requires one to associate a unique reactant conformer with each corresponding transition state conformer, MSVTST does not require this onetoone mapping; and actually based on the local equilibrium assumption in transition state theory, one should not associate a given transition state conformer with a unique path originating from a specific reactant conformer. We view this advantage of MSVTST as the correct formulation of transition state theory with multiple conformers.
3.3 Multipath variational transition state theory
MSVTST only considers one single lowestenergy path of the reaction. The multistructural effect of the saddle point (i.e., the multiple conformers of the saddle point) can lead to different reaction paths for the same chemical reaction (with different barrier heights). A more complete approach would be to explicitly consider the contributions from all these different reaction paths (of the same reaction) to the total reaction rate constant, and treat the variational effects and quantum tunneling of each reaction path individually. Multipath variational transition state theory^{149,150} (MPVTST) is a generalization based on MSVTST.
We now derive the MPVTST rate constant for a unimolecular reaction (for a bimolecular reaction the derivation follows a similar manner). For a unimolecular reaction, the MPCVT rate constant is the sum of the rate constants of all the reaction paths, which includes all the contributions from different paths into the total reactive flux:

 (129) 
where
V^{‡}_{1} is the classical barrier height of the lowestenergy reaction path;
Q^{‡SST}_{p} is the singlestructural rovibrational torsional partition function for the saddle point of path
p, which represents the contribution from the
pth saddle point structure to the total conformational rovibrational torsional partition function of the transition state:

 (130) 
where
U_{p} is the relative electronic energy of the saddle point
p with respect to the zero of energy set to be the electronic energy of the lowestenergy saddle point structure. Notice that the differences between the barrier heights of various paths have already been included in the Boltzmann factor of the singlestructural rovibrational torsional partition function of the saddle point of path
p in
eqn (130), and therefore in
eqn (129) only the
V^{‡}_{1} of the lowestenergy path was explicitly written out.
Eqn (129) can be further modified as follows

 (131) 
and it is equal to

 (132) 
where
Q^{‡SSHO}_{1} is the singlestructural rovibrational partition function of the saddle point of the lowestenergy path. Notice that if we included all the reaction paths, then the denominator of the first term in
eqn (132) is equal to the MST conformational rovibrational partition function
Q^{‡MST}. Finally,
eqn (129) can be expressed in the following form

k^{MPCVT/T}(T) = 〈γ〉F^{MST}_{fwd}k^{TST}_{1}  (133) 
where
k^{TST}_{1} is the conventional (singlestructural) TST rate constant computed based on the lowestenergy path. The averaged generalized transmission coefficient 〈
γ〉 is defined as

 (134) 
The multistructural and torsional anharmonicity factor
F^{MST}_{fwd} is

 (135) 
As we can see from
eqn (133), if we only include the single lowestenergy path, the MPVTST rate constant reduces to the MSVTST rate constant.
In the above approach, the torsional anharmonicity is only included for the saddle point and the reactant(s); the torsional anharmonicity at the variational transition states are approximated as its values at the saddle point. We can also treat the torsional anharmonicity as a function of reaction coordinate s along the MEP by replacing the singlestructural HO rovibrational partition functions Q^{XSSHO} (X is saddle point, reactant, or variational transition state) in the above equations with the singlestructural rovibrational partition function with torsional anharmonicity Q^{XSST};^{151} the computation of variational transmission coefficient is also based on Q^{XSST}. For the generalized transition state, Q^{GTSST}(s) treats the torsional anharmonicity as a function of the reaction coordinate s. The location of the dividing surface (or equivalently speaking, the value of variational transmission coefficient Γ^{CVT}_{p}) will be different from the previous approach, because the generalizedtransitionstate Gibbs free energy of activation ΔG^{GT,0}_{act}(T,s) is now evaluated with the torsional anharmonicity; and hence the position of the generalizedtransitionstate dividing surface, where ΔG^{GT,0}_{act}(T,s) reaches the maximum, could be different, i.e., the locations of the variational transition states are possibly different from the ones that are determined based on Q^{XSSHO}. As a further complete treatment, in the full MPVTST, we use Q^{XMST}(s) to compute the variational transmission coefficient, which is computationally prohibitive.
The microcanonical variational transition state theory can also be generalized to include MST effects; this is multistructural microcanonical VTST (MSμVT). And this is done by

 (136) 
in which the number of states for the saddle point is computed from the integration of density of states, and the density of states is computed from the inverse Laplace transform of the
Q^{MST}.
3.4 Truncation of paths included in MPVTST
Ideally, all the reaction paths generated by the multiple conformational structures of the transition state structure should be included in the computation of MPVTST. However, for even a mediumsized chain molecule, the number of conformers of the transition state structure is very large, and it is simply unaffordable to perform VTST calculations on all of the paths. In eqn (134), formally, the summation is from p = 1 to the total number of distinguishable reaction paths P (including paths that are mirror images), which is equal to the number N of distinguishable structures of the transition state (TS). Practically, we have to limit the number of paths included in MPVTST calculation, that is to say, in the practical MPVTST calculation, P < N. If P = N, then the denominator of eqn (134) is exactly equal to Q^{MST}(TS); and if P < N, then we have introduced truncation errors in the computed MPVTST rate constant.
A recent study^{152} compares the MPVTST results including different numbers of paths to the MPVTST with all the paths included; the comparison is for the hydrogen abstraction reaction of tertbutanol. In total there are 46 reaction paths. It was found that the variational transmission coefficient could be very different from path to path; for the lowest energy path, the variational transmission coefficients at all temperatures are close to unity, while for high energy paths, the variational transmission coefficient decreases almost linearly as temperature increases. The tunneling transmission coefficients are also pathdependent, and higherenergy paths usually have larger tunneling transmission coefficients. A simple hypothesis to explain this trend is that taller barrier tend to be thinner barriers; interestingly, this explanation does not support the observations in this study, where we examines the widths of the V^{G}_{a} potential curves for different paths in the energy range that has significant tunneling contributions, and we found that the barriers have very similar widths. Therefore we must attribute the difference to something else, and we concluded that because the higherenergy paths can tunnel at energies farther below the barrier top, they are able to gain more from tunneling contributions. As for the error introduced by truncating the number of paths, it was found that at lower temperatures (less than 300 K), one typically needs to include more paths (about 30% of the total number of paths) in order to have an error smaller than 20%, whereas at higher temperatures, including only a few paths in MPVTST is good enough (as compared to MPVTST with all paths being included).
4. Loose transition states
We can divide transition states into two classes. In one class, the transition state is loose. A loose transition state is composed of fragments that rotate with respect to one another freely or almost freely in more than one dimension (note that the existence of torsions, which are onedimensional internal rotations, is not enough to make a transition state loose). Typical loose transition states are found in bond fission reactions that do not have an intrinsic barrier, i.e., where the potential energy is in monotonically uphill or monotonically downhill in the reverse association.
In the second class, the transition state is tight, which means it is not loose as defined above. Reactions with barriers (saddle points) usually have tight transition states. Furthermore, if a saddle point is present but it is too low in energy (e.g., below the reactants of the downhill direction), then the dynamical bottleneck might be a tight transition state near the saddle point, but it also might be a loose transition state leading into a well between the reactants and the saddle point. In fact, there are probably two dynamical bottlenecks. (One can also find successive transition states, i.e., transition states in series, when there is no barrier between reactants and products.^{153}) If one of them has a much higher free energy of activation than the other, then that one should be used for transition state theory. If the two free energies of activation are comparable, one can use the unified statistical theory or the canonical unified statistical theory, as discussed in Section 4.2.
The implementation of VTST that we have reviewed so far is for tight transition states, and this version of VTST is called reactionpath VTST (RPVTST) because we search along a reaction path for variationally optimal dividing surfaces that are transverse to a reaction path; the accessible coordinate space can be described in terms of nearly harmonic vibrations and torsions around the reaction path, and the reaction coordinate is usually the distance along a reaction path, typically the MEP. Reactionpath VTST is not usually valid for a loose transition state because the accessible coordinate has wideamplitude motion that cannot be well descried by torsions and nearly harmonic stretches and bends. It is still possible to use reactionpath VTST for barrierless reactions by computing a reaction path without starting from the saddle point; however, this strategy is expected to give realistic results only (if at all) at fairly high temperatures where the transition state may tighten up. A more appropriate way to treat loose transition states is with loose transition state theory,^{154} which has now been well developed in a form called variablereactioncoordinate VTST^{155–159} (VRCVTST), as discussed next.
4.1 Variablereactioncoordinate variational transition state theory
Next we review the VRCVTST theory. The VRCVTST methodology is available in the Polyrate code^{42} and also in the earlier Variflex code.^{160} The codes were written completely independently. The present article discusses the way that the calculations are done in Polyrate.
In VRCVTST, the vibrational modes are classified as conserved modes and transitional modes. The frequencies of the vibrational modes of reactants are assumed to be conserved in the association reaction. The transitional modes consist of all of the vibrational modes except for conserved modes and the overall rotational and translational modes. During the association, the transitional modes are converted from the free rotational modes and translational modes of the infinitely separated reactants to vibrational modes and overall rotation of the associated species.
The highpressurelimit rate constant given by VRCVTST for association (R_{1} + R_{2} → P) is

 (137) 
where
g_{e} is the ratio of the electronic partition function of transition state to the product of the electronic partition functions of reactants; the ratio of the rotational symmetry numbers
σ_{1}σ_{2}/
σ^{‡} is equal to 0.5 if the reactants are the same, and it is equal to 1 if the reactants are different;
Q_{1} and
Q_{2} are the rotational partition functions excluding rotational symmetry number for the reactants;
μ is the reduced mass for the relative translational motion of the associating reactants;
N(
E,
J,
s) is the number of accessible states at energy
E and total angular momentum
J, at reaction coordinate
s.
The meaning of the reaction coordinate s in VRCVTST is different from the one in reactionpath VTST. In VRCVTST, the reaction coordinate is defined by pivot points. Pivot points are used to orientate the two fragments in the VRCVTST calculation. The locations of the pivot points are critical to the final computed rate constant; they should be chosen to yield the variationally lowest rate constants at each of the temperatures that one considers, and this is generally done by trial and error. Fixing the location of pivot point at the center of mass of the fragment or at the atomic nucleus of the dissociating/forming bond, does not necessarily lead to variationally lowest rate constant. Optimal locations of pivot points can be temperaturedependent. The reader is referred to the SI of our paper^{161} on 2CF_{2} ⇌ C_{2}F_{4} for details of the way we use pivot points to define a reaction coordinate.
In the single faceted VRCVTST, each reactant has one single pivot point, and we denote their coordinates as P_{1} and P_{2}; the reaction coordinate s is defined as s = P_{1} − P_{2}. Depends on the location of the pivot points, s may or may not be (and generally it is not) the distance between the two centers of mass of the two reactants. In multifaceted VRCVTST, each reactant can have more than one pivot points. We denote the distance between one pivot point on reactant 1 and another pivot point on reactant 2 as r_{ij}, where i is the index of pivot points on reactant 1, and j is the index of pivot points on reactant 2. The reaction coordinate s is defined as the minimal value of the r_{ij} set:
According to variational transition state theory, the pivot points should be defined so that the integral in eqn (137) is minimized.
The number of states N(E,J,s) is computed as

N(E,J,s) = 〈N_{q}(E,J,s)〉_{Ω}  (139) 
where the bracket means the average over all possible orientations of the two reactants and of the orientation of the vector
s. We define the vector
s as a vector connecting two pivot points, one in each fragment. The subscript
q means a specific orientation of the whole system. And this average is carried out by Monte Carlo sampling, where the orientations of the system are randomly sampled over the whole phase space corresponding to a given
s. The practical implementation of the algorithm is presented in detail in the SI of our paper on 2CF
_{2} ⇌ C
_{2}F
_{4}.
^{161} The number of states at a given orientation is computed as

 (140) 

 (141) 
where the bracket denotes the averaging over all possible orientations of the vector of total angular momentum
J while its length is fixed.
v =
v_{1} +
v_{2} + 3, where
v_{1} and
v_{2} are the number of orientations sampled for reactant 1 and 2;
I^{(k)}_{i} is the
ith moment of inertia of the
kth (
k = 1, 2) reactant, and
I_{i,q} is the
ith moment of inertia for the system as a whole at a given orientation
q.
n^{(12)} is the unit vector pointed from one pivot point on reactant 2 to the one on reactant 1;
d^{(k)} is the vector connecting the center of mass of the
kth reactant to its pivot point;
n^{(k)}_{i} is the unit vector directed along the
ith principal axis of the
kth reactant; and
k_{a} is number of monoatomic reactant, which is equal to 0 or 1.
As in reaction path VTST, the variational calculation can be done in different ensembles. For VRCVTST, the optimization of the transitionstate dividing surface can be done by canonical variational transition state theory (CVT), microcanonical variational transition state theory (μVT), or energy and total angular momentum resolved microcanonical variational theory (E, JμVT). In CVTVRCVTST, N(E,J,s) is calculated for all E and J at a given value of s, then the rate constant k(T,s) is minimized with respect to s at each T. In microcanonical μVTVRCVTST, the integral over J in eqn (137) is computed first to obtain N(E,s); and this is done by integrating eqn (140) over J and then averaging N_{q}(E,s) viaeqn (139). Then N(E,s) is minimized with respect to s at each E, and this optimized N(E) is used for the integral over E to obtain the final rate constant. In E, JμVT, N(E,J) is minimized with respect to s for each E and J, and then the integral is carried out to obtain final rate constant; E, JμVT yields the smallest rate constant as compared to the former two.
4.2 Canonical unified statistical theory and related extensions
Here we consider the case of two or more transition states in series. This should not be confused with multiple transition states in parallel, which are treated by multistructural variational transition state theory, multipath variational transition state theory, or ensembleaveraged variational transition state theory.
We will limit our discussion here to two transition states in series, but the generalization to three or more is straightforward. When there are two transition state regions in series, we may call them inner and outer transition state regions.^{162} The outer region is where the centrifugal barrier is located, and in some cases it may be treated using only the longrange form of the interaction potential, although in general – using the language of association reactions – one must also consider the transformation of longrange rotations to librations (which is a name for loose vibrations) and eventually to tight vibrations as the fragments approach. The inner region is where chemical interactions (valence forces) have overtaken the longrange forces. In some cases, the inner region provides the highest free energy barrier at high temperatures, and the outer region does so at low temperatures.
The canonical unified statistical theory^{22,163} (CUS) has been developed as an extension to canonical ensembles of Miller's unified statistical theory^{164,246} (US) for microcanonical ensembles, which is based on the branching analysis of Hirschfelder and Wigner.^{245} In this theory, local fluxes are calculated for dividing surfaces at both bottlenecks of the chemical process and also at a dividing surface between them. This provides the simplest way to allow for a statistical branching probability for the reaction complex to redissociate to reactants. In the case of two transition states in series, the CUS rate constant is:

k^{CUS}(T) = k^{CVT}(T)Γ^{CUS}(T)  (142) 
where
Γ^{CUS}(
T) is the CUS recrossing transmission coefficient, which plays the role of including additional recrossing corrections on the CVT rate constants. The
Γ^{CUS}(
T) factor is calculated as:

 (143) 
where
q^{CVT}(
T) is the rovibrational partition function for the canonical variational transition state (the location of which is denoted as
, which corresponds to the highest maximum point on the generalizedtransitionstate Gibbs free energy surface),
q^{max}(
T) is for the second highest maximum point on the generalizedtransitionstate Gibbs free energy surface (the location of which is denoted as
), and
q^{min} is for the lowest local minimum point that is between
and
on the generalizedtransitionstate Gibbs free energy surface. All these three partition functions are computed based on the zeroofenergy being at the minimum of reactants’ potential well (
i.e., the value of
V_{MEP} at
s = −∞); notice that, there is a relation between the generalizedtransitionstate vibrational partition function
Q^{GT}(
T,
s) that is based on the zero of energy being at the value of
V_{MEP}(
s) of the potential on MEP, and the vibrational partition function
q^{GT}(
T,
s) that is based on the former definition of zero of energy:

q^{GT}(T,s) = exp[−V_{MEP}(s)/k_{B}T]Q^{GT}(T,s)  (144) 
Clearly, if , , and are the same, as in the case which there is only one local maximum on the generalizedtransitionstate Gibbs free energy surface along the reaction coordinate, the CUS recrossing transmission coefficient Γ^{CUS}(T) = 1.
The CUS theory has been tested for various systems against accurate quantum dynamics rate constants or the best available experimentally measured rate constants.^{165–170,361} It does not necessarily get the recrossing right, but it does provide a smooth connection between the regime where the inner bottleneck is dominant and the one where the outer bottleneck is dominant.
The CUS approximation is valid at high enough pressure where the statistical dissociation of the complex is promoted by collisioninduced energy redistribution, and it is valid at low pressure if intramolecular vibrational relaxation is fast enough. But when tunneling is important, it has the same issues with total energy conservation as those discussed in relation to CVT in Section 2.3.4.
The CUS rate constant k^{CUS} can also be expressed by using the rate constant (or the reactive flux coefficient) of passing the highest variational transition state k^{I}, the second highest variational transition state k^{II}, and the formation of the lowest local minimum complex k^{C}, on the Gibbs free energy surface along the MEP, as follows:

 (145) 
And when there exists multiple parallel branching reactions after the complex, then it is straightforward to generalize eqn (145) to the following:

 (146) 
where
k^{II}_{i} is the rate constant for the
ith parallel reaction after the complex; here the transition state theory gives the overall reactive flux. And this is called the competitive canonical unified statistical (CCUS) model.
^{248} The CCUS model has been widely used in studying reaction process with competitive reaction paths
^{171–177,248} with the VTST theoretical framework, that is with the statistical model. One especially interesting application is the use of CCUS to explain the experimental KIEs for the S
_{N}2 ion–molecule reactions,
^{247} and for the competition between E2 and S
_{N}2 reactions.
^{248} A nonstatistical extension of CCUS model will be discussed in Section 6.
5. Pressuredependent reaction rate constants
5.1 Overview
For reactions with two reactants and more than one product, it is usually assumed that collisional repopulation of states depleted by reaction is fast enough to maintain the thermal equilibrium distribution of states, and hence the rate constants for bimolecular reaction can be treated as pressure independent. However, in general we must consider pressure dependence when either the reactant or the product (or both) is unimolecular, i.e., association reactions, unimolecular dissociation, and unimolecular isomerization. Variational transition state theory provides the highpressurelimit rate constants. However, in practical laboratory measurements, the pressure is seldom high enough to measure the highpressure limit, and measured rate constants can be significantly different from the ones computed for the highpressurelimit. For thermal unimolecular reactions, the rate falls off as pressure decreases, and this is called the falloff effect. It is a consequence of a nonequilibrium state distribution in the reactant.
A master equation^{34,178–182} may be used to accurately treat nonequilibrium effects if the full matrix of stateselected reaction rates and statetostate energy transfer rates is available.^{34,181,182} This is essentially never the case, but a master question can still be used if one makes additional assumptions. For example, one may assume that the unimolecular statespecific rate constants do not depend on the full state specification but rather only on the total energy or on the total energy E and total angular momentum J. We will make the former assumption. Then the master equation becomes onedimensional, and it may be written:^{183,184}

 (147) 

 (148) 
where
p_{i}(
E,
t) is the population of isomer
i with energy
E at time
t,
ω_{i}(
T,
P) is the collision frequency (which is in the unit of s
^{−1}; it is the bimolecular collision rate constant times the concentration of bath gas M),
P(
E,
E′,
t) is the probability density of collisional energy transfer from energy
E′ to
E,
y_{nA}(
t) and
y_{nB}(
t) are concentrations of the
nth arrangementchannel reactants with A and B distinguishing between the two reactants,
k_{ij}(
E),
k_{in}(
E), and
k_{ni}(
E) are the microcanonical rate constants for isomerization, association and dissociation reactions,
b_{n}(
E,
T) is the Boltzmann distribution for bimolecular reactant channel
n, and
N_{isom},
N_{react}, and
N_{prod} are the numbers of isomers, bimolecular reactant channels, and bimolecular product channels, respectively.
Popular master equation solvers are available, such as Mesmer,^{185}Multiwell,^{186} and RMG.^{187} The computations of microcanonical rate constants in these codes do not explicitly include variational effects, microcanonical multidimensional tunneling, or multistructural and coupled torsional potential anharmonicity However, in some cases variational effects are included effectively, for example, many Mesmer calculations have been performed for loose TSs with good success by first computing the parameters for a generalized threeparameter Arrhenius expression by using Variflex, which includes variational effects, and then using the inverse Laplace transform method for the generalized 3parameter Arrhenius expression to compute k(E). The inverse Laplace transform method can be viewed as an approximate way to include various dynamics effects if the computed canonical rate constants include those effects. Multiwell calculations have also been reported using variational effects; Multiwell includes a code (Ktools) for explicitly calculating k(E,J) (the rate constants as a function of energy and total angular momentum) for loose TSs (even those with more than one bottleneck); Ktools also computes canonical VTST rate constants in two ways: based on the k(E,J)s and based on a canonical VTST approach (however, in Multiwell, the canonical variation is carried out based on computing a large number of trial TST rate constants at a given temperature along the reaction coordinate, and its generalized normal mode analysis is based on Cartesian rather than curvilinear coordinates, which, as we discussed in Section 2.5, are not well suited for VTST and should be used with caution). Neither Mesmer nor Multiwell has the capability to automatically compute multistructural effects and the effects of coupled internal rotors.
However, even when adequate dynamical methods can be employed, master equation computations add an extra layer of complexity, and it would require a significant effort to compute multidimensional tunneling by summation over calculations for a large number of vibrationally excited states for each local conformational structure. In many cases we would like to use approximate chemical models to predict the pressure dependence without solving a master equation.
Pressuredependent elementary reactions can be classified, based on their mechanism, as chemically activated and thermally activated. In a chemical activation mechanism, a reactive adduct is generated via a chemical reaction; and in a thermal activation mechanism, the rovibrationally excited unimolecular states are generated via thermal collisions between the reactant and a bath gas. In the chemical activation mechanism, because collisional energy transfer to stabilize the energized adduct is slow at low pressure, the lowpressure rate constants of the formation of the stabilized adduct are lower than the highpressurelimit, and the rate constants of the further isomerization/dissociation of the formed adduct are larger than the highpressure equilibrium rate constant. Experimentally, Rabinovitch and coworkers are among the pioneers who carried out various measurements in investigating the kinetics of chemically activated reaction systems.^{188} However, in thermally activated unimolecular reactions, the lowpressure rate constants are smaller than the highpressure equilibrium rate constant because collisional energy transfer is too slow to repopulate the reacting states of the reactant (i.e., to form the energized unimolecular reactive states), and hence – as already mentioned – the pressure effect for this kind of reaction is called “falloff”. In both mechanisms, the microcanonical rate constants can be obtained via quantum RRK (QRRK) theory^{189} or via the full microcanonical variational transition state theory (sometimes called variational RRKM). Barker, Westmoreland, Dean, and Bozzelli and their coworkers^{190–197} have carried out important pioneering works in the developments and applications of QRRK theory. In many approaches, QRRK theory has been replaced by RRKM theory (i.e., transition state theory) or full variational transition state theory, which are more complete theories, but QRRK theory can significantly lower the computational effort as compared to full microcanonical VTST with multidimensional tunneling and anharmonicity. The recently proposed systemspecific quantum RRK theory^{78,198,199} (SSQRRK) is able to incorporate variational effects, multidimensional tunneling, and various anharmonicity effects into the microcanonical rate constants, without any additional cost except for the computations of the highpressurelimit rate constants. In conventional RRKM theory for k(E), energyresolved variational effects, multidimensional tunneling and multistructural and coupled torsional potential anharmonicity are not included; some of these dynamic effects can be directly incorporated into microcanonical rate constants in RRKM theory only with considerable effort. SSQRRK can be viewed as an effective way of recovering the abovementioned effects in k(E) from the highpressurelimit canonical k(T) that do include those effects in the canonical ensemble; and we use it to compute the pressuredependence so that the computed k(T,p) can connect with the correct highpressurelimit.
5.2 Thermal activation
The thermal activation mechanism for unimolecular reaction X → P is 
 (149) 

 (150) 
where X is the reactant, M is bath gas (such as He, Ar, N_{2}, etc.), P is the product, and X* is the rovibrationally excited (energized) complex, which is generated by collisions with a thermal ensemble of bath molecules M. (Of course M could be X, i.e., reaction could be occurring in neat X, and the equations can easily be modified to treat that special case, so we need not discuss it explicitly.) The parenthetical arguments following X and X* denote that X is thermalized at temperature T, and X* has internal energy E, which is greater than the threshold energy E_{0} for the reaction step. (M is also thermal at temperature T, and the translational energy of X* is thermal at temperature T, but these conditions are not explicit in the notation.) The rate constant for the thermal activation step is k_{act}, and we will call the reverse step deactivation and label it k_{deact}, although in this step the decrease of energy in X* could still yield a molecule that is activated enough to react, except in the limit of so called strong collisions, in which each collision fully deenergizes X*. The rate constant for the step in which reactant is converted to product is called k_{con,X} (this could be either isomerization or dissociation). The phenomenological unimolecular reaction rate constant is defined as 
 (151) 
In the classical Lindemann theory, all of the rate constants are considered to be energyindependent. If we use steadystate approximation to treat X*, i.e., if we assume

 (152) 
then we get

 (153) 
If we assume idealgas behavior, [M] can be replaced by p/RT, where p is pressure. Then, since k_{uni}[X] = k_{con,X}[X*], we have

 (154) 
In the highpressure limit, the Lindemann unimolecular rate constant becomes

k^{∞}_{uni} = k_{act}k_{con,X}/k_{deact}  (155) 
and the unimolecular reaction is of pseudofirst order; whereas in the limit of very low pressure, the lowpressurelimit
k^{0}_{uni} is much smaller than
k^{∞}_{uni}:

k^{0}_{uni} = (p/RT)k_{act}  (156) 
which depends linearly on the pressure, and the reaction is of second order. The falling off of the
k_{uni} rate constant predicted by classical Lindemann theory starts too early,
i.e., the computed pressure at which significant falloff effect begins is too high as compared to experiment.
For higher accuracy, the energy dependence of k_{con,X} should be considered. One could use microcanonical VTST to calculate the energydependent rate constants, just as RRKM theory uses microcanonical conventional TST energydependent rate constants. Alternatively one could calculate thermal rate constants and obtain the fixedenergy ones by an inverse Laplace transform as mentioned at the end of Section 2.1, but inverse Laplace transform is often unstable. What we will do instead is to use SSQRRK as an effective kineticsspecific inverse Laplace transform to convert canonical k_{con,X} rate constants to microcanonical ones; SSQRRK is a physical model for achieving this goal, but without explicitly carrying out the inverse Laplace transform in its actual computation.
In the thermal activation step, the thermally equilibrated reactant X at temperature T collides with bath gas M to produce X* with internal energy E, and a key assumption of classical RRK theory, which is retained in QRRK and SSQRRK, is that we neglect rotational energy and assume the vibrational energy E is randomly distributed among all the modes (nevertheless, in SSQRRK the rotational contribution is implicitly included in k(E) because the canonical highpressurelimit includes rotational modes, and thus they are implicitly in A^{K} and E_{0}, which are discussed below). For simplicity we treat X and X* as multimode harmonic oscillators with all frequencies the same (i.e., an effective geometric mean frequency); this will not be as serious an approximation as it might seem because the resulting formula will be used only to make a kineticsspecific inverse Laplace transform of a thermal rate constants calculated without this assumption. In calculations, we will use the geometric mean frequency of X, and we will call in (in wavenumbers). With all the internal energy distributed among harmonic oscillators of the same frequency, we may change the independent variable from internal energy E to number of vibrational quanta n:

n = E/hc  (157) 
where
c is speed of light. Following a similar derivation to that in classical Lindemann theory, but employing
ndependent rate constants, the unimolecular rate constant is

 (158) 
in which
m is the number of quanta at the threshold energy
E_{0},

m = E_{0}/hc  (159) 
and the determination of threshold energy in SSQRRK will be discussed later.
The ratio of the thermal activation and deactivation rate constants in eqn (149) is a mixedensemble equilibrium constant, representing the thermal equilibrium of species X* in a microcanonical ensemble at energy E with species X in a thermal ensemble with temperature T. In the QRRK model, this is the fraction of molecules at temperature T that have n quanta of excitation, which is given by^{200}

 (160) 
where
s is the number of vibrational degrees of freedom of X* (which is 3
N − 6 for a nonlinear molecule). We assume that
k_{deact} does not depend on
E,
i.e., does not depend on number of vibrational quanta
n, but does depend on
T; then we solve
eqn (160) for
k_{act}, which yields

 (161) 
In QRRK, a unimolecular reaction happens when a specific vibrational mode associated with the reaction coordinate possesses threshold energy E_{0}. Therefore the QRRK rate constant equals a frequency factor A^{K} (the superscript stands for Kassel) times the fraction of molecules that have at least m quanta of excitation; this is given by^{181,201}

 (162) 
Substituting eqn (159), (161) and (162) into (158) yields the QRRK rate constant, but in order to finish the calculation we must specify A^{K}, E_{0}, and k_{deact}. The computation of k_{deact} will be discussed in Section 5.4, and next we turn our attention to A^{K} and E_{0}.
Taking the highpressurelimit of eqn (158), using eqn (159), and summing the series yields the following pure Arrhenius form:

k^{∞}_{uni}(T) = A^{K}exp(−E_{0}/k_{B}T)  (163) 
and this motivates the approach we proposed to obtain the parameters needed in SSQRRK computations. The first step is to calculate the highpressurelimit canonicalensemble VTST rate constant (for example, an MPCVT/SCT calculation that includes variational effects, multidimensional tunneling, various anharmonicity). Next, we fit the canonicalensemble highpressurelimit rate constant by the fourparameter formula
^{202} 
 (164) 
where
B,
T_{0},
E′, and
n′ are fitting constants. In combustion community, the widely used fitting formulas are either conventional Arrhenius formula or the threeparameter fitting formula in the form
BT^{n}exp(−
E/
RT), which assumes a linear dependence of activation energy on temperature; however, this linear dependence can be very inaccurate at low temperatures, where the tunneling nonlinearly reduce the activation energy significantly,
^{150} and also at high temperature. Our fourparameter formulas are more physically sound, and they permits a nonzero rate constant for an exothermic reaction in the lowtemperature limit; furthermore they are more accurate in terms of the fitting error.
^{202}
Based on our fourparameter fitting formula, the local Arrhenius activation energy (which may also be called the Tolman activation energy)^{203} for temperatures near a given temperature T′ is obtained from the negative slope of the Arrhenius plot,

 (165) 
which yields

 (166) 
With this energy of activation, we have a local Arrhenius fit for temperature T given by
k^{∞}(T) = A(T)exp[−E_{a}(T)/RT] 
with
A(
T) obtained by

A(T) = k^{∞}(T)exp[E_{a}(T)/RT]  (167) 
Then E_{a}(T) and A(T) are used for E_{0} and A^{K} in the QRRK expression given by eqn (158), (159), (161) and (162). This yields SSQRRK. Notice that the parameters E_{0} and A^{K} are constants in conventional QRRK, but they depend on temperature in SSQRRK, which is a major change. This procedure allows one to incorporate the full apparatus of VTST (including variational optimization of the transition state, tunneling, and anharmonicity) in the highpressurelimit rate constants and approximately transfer it to the microcanonical energydependent rate constants, without actually doing the microcanonical computations. We view the subsequently obtained SSQRRK microcanonical rate constants as an effective approximation to the microcanonical variational transition state theory rate constants including variational optimization, tunneling, and anharmonicity.
If one wants to do the microcanonical computations (which could be expensive in order to fully account for MST anharmonicity and microcanonical multidimensional tunneling), the sum in eqn (158) can be replaced by an integration, which yields

 (168) 
in which

 (169) 
where
ρ(
E) is MST rovibrational density of states of X,
Q^{X}(
T) is the MST rovibrational partition function of reactant X,
E_{thr} is the threshold energy. Notice that, in this rigorous microcanonical treatment where the microcanonical rate constants are computed by μVT theory (which can be viewed as a variational RRKM theory), the threshold energy is the ZPEincluded electronic potentialenergy barrier height; this is different from the way of determining the threshold energy in SSQRRK theory, which is a temperaturedependent activation energy. Also, the rotational contributions are explicitly treated in this formalism. The microcanonical rate constants in
eqn (168) can be computed by MSμVT with multidimensional tunneling. In the case of smallcurvature tunneling (SCT), the MSμVT/SCT rate constant is

 (170) 
where the cumulative reaction probability
N^{‡,MST/SCT} is the convolution of the SCT tunneling probability
P^{SCT}(
E) with the density of states of the saddle point, in which the summation goes over all rovibrational states of all conformations of the transition state, from vibrational groundstate and continues up in energy until it converges or until one reaches the dissociation energy of the transition state. This is computationally demanding even for mediumsized systems.
Finally, let us consider use SSQRRK for treating unimolecular dissociation with multiple parallel dissociation reactions.

 (171) 

 (172) 
where the reactant X has multiple parallel dissociation reactions, of which the products are respectively P
_{1}, P
_{2},…, P
_{N} (notice that here the notation P
_{i} for the
ith dissociation reaction represents the product(s) of this dissociation reaction). The overall phenomenological unimolecular reaction rate constant (for the total dissociation rate constant of the reactant X) is defined as:

 (173) 
which, by using the condition [X] = [P
_{1}] + [P
_{2}] + ⋯ + [P
_{N}], can be rewritten as the summation of the rate constants of all the parallel dissociation reactions

 (174) 
where

 (175) 
And the SSQRRK pressuredependent rate constant for the ith dissociation reaction, which is a straightforward generalization of eqn (158), is:

 (176) 
where
m_{i} is the threshold number of quanta for dissociation reaction
i, which is determined by using
eqn (159), and
k^{(i)}_{con,X}(
n) is determined by
eqn (162).
Eqn (176) is equivalent to eqn (A5) of Dollet
et al.,
^{204} who applied it to study falloff in the cases such as SiH
_{4} → SiH
_{3} + H or SiH
_{2} + H
_{2}, although theirs work was based on conventional TST without including the dynamic effects we consider in SSQRRK.
The addition of rate constants for parallel reactions, as in eqn (174), is clearly valid for microcanonical rate constants and for highpressurelimit thermal rate constants for which the energy distribution is thermal. But for application to the falloff region, eqn (176) involves the additional assumption that energy transfer competes with energy relaxation in the same way for all product channels. But this is not necessarily the case in the falloff region if highenergy states are depleted more rapidly by reactions with lowerenergy thresholds than they are by reactions with highenergy thresholds. This phenomenon has been discussed in conjunction with the treatment of the unimolecular reactions of 2methylhexyl radicals.^{205} Thus eqn (176) should be used with caution in the falloff region if the parallel reactions have a large difference in threshold energies.
5.3 Chemical activation
For bimolecular association reaction or a twostep reaction with a unimolecular intermediate, the chemical activation mechanisms are:
(a) association reaction (Y + Z → YZ)

 (177) 
(b) twostep reaction (Y + Z → YZ → P)

 (178) 
where YZ* is the activated adduct generated by reaction between Y and Z, and YZ is the stabilized product (it is the intermediate in the case of twostep reaction), P is the product of further reaction of the intermediate YZ*. In chemical activation, the rate constant
k_{add} of the formation of YZ* is not pressuredependent because it is formed
via a bimolecular chemical reaction, of which the value is just the highpressurelimit temperaturedependent bimolecular reaction rate
k^{∞}_{add}(
T). The rate constant for the dissociation of YZ* back to the reactants Y and Z is
k_{rev}, which is treated as energydependent; the rate constant for the formation of the product(s) from YZ* is
k_{con,YZ}(
E).
Following a similar treatment as for thermally activated unimolecular reactions, and defining k_{stab} as the phenomenological bimolecular stabilization rate constant,

 (179) 
and
k_{rxn} as the phenomenological bimolecular rate constant for reaction

 (180) 
We applied steadystate approximation to the activated adduct YZ*. Then the pressuredependent rate constants given by SSQRRK theory are:
(a) for association reaction

 (181) 
(b) for twostep reaction

 (182) 
and

 (183) 
where
f(
E) is the fraction of energized adduct (YZ*) at energy
E, which is given by

 (184) 
and
K_{act,YZ} is given by
eqn (157) and (160), and
E_{0} is the SSQRRK threshold energy for YZ dissociating back to reactants Y and Z, and
k_{rev} and
k_{con,YZ} are computed by
eqn (162) with the appropriate
m (number of quanta at threshold energy) and
A^{K} values (for
k_{rev}, these parameters are computed from the reaction of YZ dissociating to Y and Z; and for
k_{con,YZ}, from the reaction of YZ dissociation to P).
Instead of using SSQRRK to approximate microcanonical rate constants, the above equations can be converted to the form needed to use microcanonical VTST rate constants (including microcanonical quantum tunneling):
(a) association reaction

 (185) 
(b) twostep reaction

 (186) 

 (187) 
where the energy distribution function

 (188) 
and in μVT, the threshold energy is denoted as
E_{thr}, and it is ZPEincluded electronic potentialenergy barrier height, which is different from the threshold energy used in SSQRRK theory (which is a temperaturedependent activation energy).
As for the thermal activation process of unimolecular dissociation with multiple parallel dissociation reactions, which we discussed in Section 5.2, here we consider using SSQRRK for chemical activation processes with the intermediate YZ undergoing multiple further parallel dissociation reactions. The mechanism is as follows:

 (189) 

 (190) 
The phenomenological bimolecular rate constant for the formation of the product(s) P_{i} of parallel reaction i in eqn (190) is defined as:

 (191) 
The SSQRRK rate constant for the stabilization of the intermediate YZ, which is an extension based on eqn (182), is:

 (192) 
where
E_{0} is the SSQRRK threshold energy for YZ dissociating back to reactants Y and Z. As we have discussed at the end of Section 5.2, if the depletions of the highenergy states by the lowthreshold and highthreshold reactions are comparable, then the current treatment for parallel channels provides a good approximation. And the pressuredependent rate constant
k^{(i)}_{rxn} for the formation of the product(s) P
_{i} is:

 (193) 
5.4 Collisional deactivation rate constant
In the strong collision approximation, the collision efficiency is unity, which means that, in the strong collision limit, each collision produces a Boltzmann distribution of the energized complex (i.e., of the activated adduct in a chemically activated system, and the rovibrationally excited unimolecular states in a thermally activated system), which usually has a much lower average energy than the energy prior to collision (i.e., deactivation has occurred). Note that deactivation refers to reducing the total energy below the threshold for reaction. However, in reality it may takes multiple energy transfer collisions to fully deactivate the energized complex. At lower pressure, where there is more time between collisions, the energy relaxation time may be slower than the reaction time. Troe developed approximate way to treat singlewell, singlereactionchannel systems by a masterequationbased weakcollision model,^{206–209} which is derived from the full master equation with an assumption that the collisional energy transfer probability has an exponential form; Troe's collision model has been tested against various bath gases, and its performance is well validated.^{215} In Troe's model, the bimolecular collisional rate constant k_{deact} (in units of cm^{3} molecule^{−1} s^{−1}) is computed as: 
k_{deact} = β_{c}Ω_{2,2}*k_{HS}  (194) 
where k_{HS} is the hardsphere collisional rate constant, β_{c} is the collision efficiency (dimensionless), and Ω_{2,2}* is the reduced collisional integral, for which the following approximation is available:^{206,207} 
 (195) 
where ε_{A–M} is the LennardJones interaction parameter between molecule A and M and is computed as the geometric average of ε_{A–A} and ε_{M–M}. The collision efficiency is computed as: 
 (196) 
where 〈ΔE〉 is the value of the average energy transferred during both energization and deenergization processes (this quantity is also denoted as 〈ΔE〉_{all}, in order to be distinguishable from the energy transferred only for deenergization 〈ΔE〉_{down},^{179,210,211} which is a positive number) and 
 (197) 
is the thermal population of internal energies of unimolecular states above the threshold energy of the reactant normalized by a density of states factor at the threshold energy.^{212} If only the 〈ΔE〉_{down} parameter is available, then the collision efficiency could be computed as:^{179} 
 (198) 
The F_{E0} factor can be computed directly from the density of states (which is the inverse Laplace transform of the rovibrational partition function), and in this case, the multistructural torsional anharmonicity could be included in the factor F_{E0} by using the MST partition function, which is reviewed in Section 3.1, in the inverse Laplace transform; if the single structure quasiharmonic oscillator model is a good approximation, F_{E0} can also be computed by direct counting algorithms, or (with much less computational cost) the empirical Whitten–Rabinovitch approximation,^{213} which is based on the classical harmonic oscillator model. In principle, the value of the energy transfer parameter 〈ΔE〉_{all} could be determined from theoretical trajectory calculations, as done in ref. 214–216 (although the way to extract information from trajectory calculation is not completely unambiguous). Most applications use a value that has been used in the literature for similar systems; or this quantity is simply treated as a fitting parameter to reproduce the experimental results in a certain range of temperature and pressure, and then one hopes that the soobtained pressuredependence model becomes predictive for other temperatures and pressures.
For the numerical value of F_{E0}, Troe studied a set of small molecules and concluded that their values are all close to unity.^{212} Using eqn (197), we also confirmed these observations for some cases, although, for larger sized molecules and at high temperatures, F_{E0} can deviate from 1 or 1.15 by quite a large amount; in fact, for a largesized molecule with lots of vibrational degrees of freedom at high temperatures, F_{E0} grows exponentially. For instance, for C_{2}F_{4} dissociating to two CF_{2},^{161}F_{E0} is computed to be 1.01 at 300 K, 1.12 at 400 K, and 1.39 at 1000 K; for (SiH_{3})_{2}SiHSiH^{−} isomerization to (SiH_{3})_{2}SiSiH_{3}^{−},^{199} it is 1.43 at 300 K, 1.63 at 400 K, and 6.54 at 1000 K; and for H addition on toluene,^{78} it is 1.31 at 300 K, 1.47 at 400 K, 4.02 at 1000 K. Therefore it is not justified to use a constant value close to unity, as has been done^{189,190} in the literature. When the density of states is obtained via our MST partition function (see Section 3.1) and F_{E0} is computed by eqn (197), the multiple structural and torsional anharmonicity (and also, the vibrational anharmonicity for highfrequency modes when a vibrational scale factor is applied) is built into F_{E0}, and therefore the previously proposed empirical internal rotor correction and anharmonicity correction^{212} to F_{E0} should not be used.
The collision efficiency β_{c} is smaller than 1 and it is significantly smaller than 1 at high temperatures for small molecules; furthermore β_{c} decreases with increasing temperature. The computed value of β_{c} clearly depends on the value of the energy transfer parameter 〈ΔE〉_{all}, and this is often the largest source of uncertainty in falloff computations. Typically one finds that larger molecules have higher values of the energy transfer parameter;^{217,218} and a larger value of the energy transfer parameter leads to a larger value of collision efficiency (i.e., a stronger collision). As mentioned above, a common practice in the literature is to use a previously used value of the energy transfer parameter for similar systems. In principle, 〈ΔE〉_{all} should depend on energy or temperature,^{218–221} but it is often approximated as a constant; to introduce the temperature dependence, an additional empirical parameter is then needed. Furthermore, even when an experimental value is available, the experiment may correspond to a different energy or temperature than the one that is needed for the kinetics calculation. Yet another complication that occurs in chemical activation mechanisms is that the available experiments do not usually identify the adduct. For example, in our calculations on addition of H atoms to toluene, we used an experimental value for disappearance of toluene, whereas what is really needed is an experimental value for producing the adduct of H with toluene. The reader is directed to the literature^{222–228} for examples of experimental values of 〈ΔE〉_{all}.
A number of collision efficiencies β_{c} for a variety of chemical systems that are obtained in the experimental work (in some cases, they are not obtained directly or solely based on experimental measurements; they are derived with some information from theoretical models) are collected in Tardy and Rabinovitch's review;^{229} for the collected collision efficiencies in their review, the authors emphasized that “much of the presentation is qualitative in nature” and they noted a “lack of quantification”. Here we give some examples of numerical values in some of the recent SSQRRK calculations for the collision efficiency β_{c}. For H addition on toluene (Ar as bath gas, with the energy transfer parameter 〈ΔE〉_{all} being 130 cm^{−1}),^{78} it is 0.2 at 300 K, 0.05 at 800 K, and 0.003 at 1500 K. For the SO_{2} and OH association reaction in the atmosphere (assuming N_{2} as bath gas, with the energy transfer parameter 〈ΔE〉_{all} being 74 cm^{−1}),^{232} it is 0.21 at 250 K, 0.18 at 298 K, and 0.14 at 400 K. For C_{2}F_{4} dissociating to two CF_{2} (with Ar as bath gas, with the energy transfer parameter 〈ΔE〉_{all} being 250 cm^{−1}),161 it is 0.4 at 300 K, 0.2 at 800 K, and 0.1 at 1500 K.
5.5 Fitting the pressuredependent rate constant k(T,p)
Functional forms have previously been proposed for fitting the pressuredependent rate constant with respect to both the temperature and pressure,^{230,231} but recently we have suggested a simpler formula for this purpose, and it has been shown that this fitting yields quite satisfactory results.^{232} The formula involves making three onedimensional fits as a function of temperature.
Step 1: fit the computed highpressurelimit rate constant k^{∞}(T) (for a bimolecular association reaction, the units are cm^{3} molecule^{−1} s^{−1}; and for a unimolecular reaction, the units s^{−1}), which has already been done by using eqn (164).
Step 2: fit the lowpressurelimit rate constant k_{0}(T), which is defined as the value of k(T,p)/[M] in the limit of p going to zero. The lowpressurelimit rate constant is a pseudothirdorder rate constant with unit being cm^{6} molecule^{−2} s^{−1} for a bimolecular association, and it is a pseudosecondorder rate constant with unit being cm^{3} molecule^{−1} s^{−1} for a unimolecular reaction. The fitting formula used in this step has the same functional form as the one used in step 1, i.e.,

 (199) 
where
A,
n,
E and
T_{0} are the four fitting parameters,
T is temperature in K, and
R is the ideal gas constant (in the unit of 1.9872 × 10
^{−3} kcal mol
^{−1} K
^{−1}).
Step 3: fit the transition pressures p_{1/2}(T), which have the units of bar. The transition pressure is defined as the pressure at which the rate constant is half of the highpressurelimit. The p_{1/2}(T) could be fit with many possible formulas, and we give one of the fitting formulas here:

 (200) 
where
a_{1},
a_{2},
a_{3},
l_{1},
l_{2},
T_{1}, and
T_{2} are the fitting parameters, in which
a_{1},
a_{2}, and
a_{3} are unitless, and
l_{1},
l_{2},
T_{1}, and
T_{2} are in K.
The final fitted pressuredependent rate constant k(T,p) is expressed in terms of the fitted k^{∞}(T), k_{0}(T), and p_{1/2}(T) (and thus it does not require any additional fitting) by the following interpolatory equation:

 (201) 
where

 (202) 
and

 (203) 
in which the Boltzmann constant
k_{B} in
eqn (202) is 1.38 × 10
^{−22} bar cm
^{3} molecule
^{−1} K
^{−1},
p is pressure in bar,
B is a derived parameter in units of bar
^{−1}, and
d is in bar
^{2}.
Eqn (201) has the following three properties, which ensure the fitted rate constants have the right highpressurelimit, lowpressurelimit, and the transition pressures:

k(T,p = ∞) = k^{∞}(T)  (204) 

k(T,p = 0) = k_{0}(T)[M]  (205) 
and

k[T,p = p_{1/2}(T)] = k^{∞}(T)/2  (206) 
The pressuredependent activation energy E_{a}(T,p) at a given pressure p = p_{0} can then be determined by numerical differentiation by employing the following formula:

 (207) 
In the highpressure limit, the above equation would yield results that are identical to the ones given by eqn (166) at T = T′.
6. Applicability of transition state theory
In this section, we discuss the applicability of transition state theory. When we say transition state theory in such a context, we mean variational transition state theory.
Classical transition state theory without a transmission coefficient involves two main assumptions.^{19,41} The first assumption of transition state theory is that the states (states in a quantal treatment, phase points in a classical one) of the reactants are in equilibrium (thermal equilibrium when we discuss temperaturedependent rate constants as in a canonical ensemble; microcanonical equilibrium when we discuss energydependent rate constants). The question arises of whether it is an additional assumption to treat the transition states by equilibrium statistical mechanics. This question has been answered, and in classical mechanical language we have the theorem:^{43,44} “When products are absent the distribution in the transition region is identical to an equilibrium distribution except that states [phase points] lying on trajectories originating in the products are missing.” Thus it is not a separate assumption to assume that states of the transition state have a Boltzmann distribution.
The second assumption of transition state theory is that the states crossing the transition state dividing surface in the direction toward the products are on trajectories that originated from the reactant region and that they will proceed to products without recrossing the transition state.^{20}
These two basic assumptions are called the localequilibrium or quasiequilibrium assumption and the norecrossing assumption. We have stated the assumptions in classical mechanical language, but they can be rephrased quantum mechanically by using Heisenbergtransformed projection operators in the flux calculations instead of the classical concept of where trajectories initiated.^{19,233} However, one cannot fully translate the classical justification of transition state theory into quantum mechanics^{234} (for example, one no longer has a practical variational bound on the rate constant) because the uncertainty relations prevent us from simultaneously fixing the reaction coordinate at a definite value and specifying the sign of the momentum in that coordinate. Thus, in order to pass from classical transition state theory to quasiclassical transition state theory, we make another assumption, and the assumption we make is that we can replace the classical partition functions by quantum mechanical ones in calculating the flux. We may call this the quasiclassical postulate since this is the step that converts classical transition state theory to what we have called quasiclassical transition state theory. (In early papers we called it hybrid transition state theory, but now we recommend the quasiclassical notation. Note also that this is called semiclassical transition state theory in the kineticisotopeeffects literature,^{235–237} but we prefer to save the word semiclassical for WKBlike treatments of tunneling. Thus when we say “semiclassical,” we refer to including tunneling, but when they say “semiclassical” in the kineticisotopeeffects literature, they mean without tunneling, which is our quasiclassical.)
Although the quasiclassical postulate is far from rigorous, it can be justified in various ways. Wigner was the first to do this; he showed that such a substitution is justified to the lowest order in ħ.^{83,238}
Modern quantum dynamics provides an even better justification for the quasiclassical postulate. A key aspect of the quasiclassical postulate is that systems cannot pass through the transition state with less than zero point energy except by tunneling. If we take the quasiclassical postulate literally, it implies that systems actually pass through the transition state not just in the quantized ground state (i.e., with zero point energy) but in all of the quantized vibrational levels (excited levels as well as the zero point level). Because the transition state is metastable (i.e., unlike stationary states in quantum mechanics, it has a finite lifetime), the quantized levels are broadened, i.e., they are quantum mechanical resonance states with finite widths related to their lifetimes. Accurate quantum mechanical calculations have confirmed that this is actually the case.^{7–10} The levels can be observed in the accurate quantum mechanical cumulative reaction probabilities; the levels have widths related to the lifetimes of individual quantized transition states. The widths can also be related to tunneling, i.e., tunneling through the barrier in a vibrationally adiabatic potential energy curve corresponds to reacting at an energy lower than the center of the broadened vibrational level, and thus it accounts for the broadening. It is beyond the scope of the present paper to go into this in more detail, but reviews are available.^{8,9}
In addition to the three basic assumptions of quasiclassical transition state theory (quasiequilibrium, no recrossing, and the quasiclassical postulate), there are additional approximations in the treatment of the transmission coefficient unless one uses eqn (98), but eqn (98) is not generally practical.
Confirmation of the quantitative accuracy of transition state theory including tunneling transmission coefficients is complicated by the uncertainties in most available potential energy surfaces. Thus there are errors in the potential energy surface and errors in the dynamics for the given surface, and these two kinds of error can cancel or reinforce. One of the few systems for which the potential energy surface is well enough known so as not to complicate tests of the dynamics is the H + H_{2} reaction, and transition state theory with tunneling transmission coefficients that include corner cutting agree well with experiment for that case.^{106,239} It is encouraging that as the potential energy surfaces have improved in accuracy in recent years, there are many other cases where VTST agrees well with experiment. Nevertheless the best tests of the dynamics are still the tests against accurate quantum dynamics when the same potential energy surface is used because then one can compare approximate and accurate dynamics on the same potential energy surface, and the difference cannot be due to inaccuracy of the surface but must instead reflect a test of the theory used to calculate the approximate rate constants. Transition state theory has fared very well in such tests. For example, in a review of 311 such rate constant comparisons for collinear and threedimensional atom–diatom reactions in the temperature interval 200–2400 K, CVT/μOMT rate constants were found to agree with accurate quantum mechanical rate constants within 22% on average.^{6} Accurate quantal rate constants are also available for H + CH_{4} → H_{2} + CH_{3} in the temperature range 200–1000 K; in this case the average deviation is only 14%.^{240}
Transition state theory is a dynamical theory because it calculates the flux through a dividing surface. It is also a statistical theory in that it provides a statistically averaged highpressurelimit rate constant of an ensemble, where the highpressure limit is invoked to main the local equilibrium of reactant states. Strictly speaking, nonequilibrium effects are not included, although sometimes one extends the theory by including a nonequilibrium factor in the transmission coefficient and still calls the result transition theory, so one must be careful to understand what an author means. When we calculate pressuredependent reactions, we still make the equilibrium assumption for microcanonical ensembles, but we no longer assume that these microcanonical ensembles are weighted the same as in a canonical ensemble (as we do for thermal equilibrium); so the theory is still applied to equilibrium ensembles but they are not the same ones as for thermal rate constants (i.e., as for canonical rate constants). RRKM theory is transition state theory applied to microcanonical ensembles for unimolecular reactions, so discussions about whether microcanonical ensembles are equilibrated during reactions apply to both transition state theory and RRKM theory.
The localequilibrium assumption is usually expected to hold well for bimolecular reactions because the translational energy, rotational energy, and vibrational energy of reactants are assumed to equilibrate on a faster time scale than the time scale for reaction. This is usually true, but there may be deviations for very fast reactants. Thus bimolecular reactions without an intermediate are usually treated as pressure independent. The participation of an intermediate in a lowtemperature reaction involving tunneling, as discussed in Section 2.3.4, may, however, cause pressure effects, and pressure effects are prominent in the cases discussed in Section 5, namely (i) chemical activation, which is a two step mechanism with the first step bimolecular and second being unimolecular and (ii) thermally activated unimolecular reactions. Section 5 discusses the reaction of chemically activated intermediates and thermally activated intermediates as a function of total energy E (one can similarly treat their rates as functions of E and total angular momentum J). The rate at which the distribution of E becomes thermalized is not a question that transition state theory answers or makes assumptions about; but the transition state quasiequilibrium assumption enters when we assume that for those species with a given E or E and J, phase space is statistically occupied subject only to these constraints. This follows if we assume fast intramolecular vibrationalenergy relaxation (“fast IVR”), and fast IVR may be considered part of the transition state theory quasiequilibrium assumption in this context. Thus transition state theory, and hence also RRKM theory, is not applicable to very shortlived complexes that do not have time for full intramolecular energy redistribution^{241} or during the induction period of a longlived complex, i.e., during that period before a steadystate phenomenological rate constant is established.
Next we turn to the norecrossing assumption.
Hase and coworkers discussed the possible breakdown of TST based on trajectory computations for ion–molecule reactions.^{242,243} As Hase points out,^{244} “the RRKM model for the complexes may be invalid” and “recrossing of the central barrier may also be important, violating the fundamental assumption of transition state theory.” One issue is that classical trajectory computations could overestimate the recrossing because classical trajectory computations do not enforce zero point vibrational energy during the trajectory. Although conservation of zero point energy is not required during dynamics, we know – as discussed above – that cumulative reaction probabilities (and hence threshold energies) are consistent with its approximate conservation. Therefore, it is not clear to what extent these recrossing effects observed in trajectory computations are due to the classical nature of the calculations. It has also been pointed though that the quantitative validity of transitionstate theory may be quite different in classical and quantum mechanics.^{26} Interestingly, a canonical unified statistical theory^{163} (which is an extension of the earlier unified statistical theory^{245,246}) has been applied to study similar problems,^{247,248} and the computed kinetic isotopic effects (KIEs) agree reasonably well with experimental observations for gasphase S_{N}2 reactions.
Consider, as an example, association reactions involving small molecules, where the system crosses a loose bottleneck. If there are not enough lowfrequency degrees of freedom to redistribute energy, the system can hit the repulsive wall behind the well, and bounce right back, and thus the capture probability is small. This fact has been known for a long time.^{26,249} For many years, there have been studies of whether energy is randomized in unimolecular decay, as assumed in RRK, RRKM, and TST theories. This is actually the central issue about the validity of RRKM theory since RRKM theory is TST with the RRK assumption that reactants are equilibrated, which is also a fundamental assumption of TST. In some cases,^{250} TST describes smallmolecule association reasonably well.
Transition state theory does not predict product energy distributions,^{251} and it predicts chemical product distributions only when they each have their own variational transition state dividing surface. Therefore it does not predict relative reaction rates that are determined by bifurcation of a reaction path after a single dynamical bottleneck that applies to both reactions.^{252–256} Similarly, if there is a highfreeenergy bottleneck followed by two parallel lower ones, and the later dynamical bottlenecks are not separated from the earlier one by a thermalized intermediate, one should not try to predict relative reaction rates from the two lower ones by transition state theory – not because transition state theory is in error but because it does not address this issue. When there is a bifurcation after the overall bottleneck, TST only gives the sum of reaction rates along the two paths. One can however, try to model such problems with beyondtransitionstatetheory assumptions. For example, we combined VTST with nonstatistical assumptions to study bifurcation in the BH_{3} + propene system, and we found that our computed selectivity agrees well with experimental observations.^{257} This method is called the canonical competitive nonstatistical (CCNM) model.
The dynamical problem that is modeled by CCNM is that of a reaction with an overall dynamical bottleneck (the maximum of the generalized transition state free energy along the MEP) that is followed by a branching into two possible products. In such models, we have two types of reaction mechanisms.^{257} The indirect mechanism, which is computed by VTST/MT, applies to some fraction of the molecules in the ensemble that become equilibrated in an intermediate complex (which is denoted as “C”), and the reaction of the equilibrated intermediate is controlled by the second set of the dynamical bottlenecks (that is the two transition states, which are denoted by 2A and 2B, leading to the two possible products); and the direct mechanism, which is modeled by nonstatistical phase space theory,^{258,259} applies to the remaining portion of the ensemble members, which keep a significant amount of energy in the reaction coordinate degree of freedom until after passing the second bottleneck region (sometimes this is called a “hot” molecule reaction). The total rate constant k_{tot}, which is the sum of the direct (k_{dir}) and indirect (k_{ind}) components, is computed as:

 (208) 
where
k_{*2} = min(
k_{C},
k_{*2A} +
k_{*2B}), and
k_{i} is the reactive flux coefficient for a dividing surface at location
i (
i = *1, C, 2A, *2B; “*1” is for the first highest local maximum, “C” is at the location of the complex, “*2A” and “*2B” are for the two branching reaction A and B). The direct component (
k_{dir}) is computed as:

k_{dir} = k_{*1}k_{*2}/k_{C}  (209) 
and the total reaction rate constants for the two branches A and B are:

 (210) 

 (211) 
where
R_{A/B} is the branching ratio of the direct reactive components, which is approximated as the equilibrium constant
K_{A/B} given by statistical phase space theory
^{260,261} with additional nonstatistical corrections:

 (212) 
where the rovibrational partition function
Q_{X} (X = A or B, which denote the final product of reaction A or reaction B) is

 (213) 
in which
Q_{X,⊥} is the rovibrational partition function that includes all the boundmode vibrations that are orthogonal to the reaction coordinate
s, but excluding the reaction coordinate
s itself, which becomes a boundmode vibration in the product (and it is denoted as mode
F, and its corresponding vibrational quantum number is
n_{F} and vibrational energy is
ε_{nF}). The parameter
τ_{relax} is the characteristic relaxation time for intramolecular vibrational relaxation, which cannot be obtained with a statistical theory, and it is treated as an empirical parameter in the CCNM model (although, some representative values, which are about 100 fs, are available
^{258} in the literature). The
τ_{nF} characterizes the duration of the collisional thermalization for a given vibrational state, and it is estimated in the model as the ratio of a characteristic distance (
s_{product}–
s_{TS1}) in coordinates scaled to a reduced mass to a characteristic velocity in the reaction coordinate for reaction coordinate quantum number
n_{F}.
^{257}
Although, we do not expect this proposed CCNM model to be able to quantitatively treat any bifurcation systems without any further modification, our point is that the VTST theory itself (when applied with a highlevel potential energy surface and appropriate treatment of anharmonicity and tunneling) can be accurate for the total rate constant but it should not be applied to predict what happens after the system passes the overall dynamical bottleneck.
Glowacki and coworkers^{262} presented an interesting study for the same BH_{3} + propene system by an alternative method. They use a master equation approach with the microcanonical rate constants needed in the model being computed by RRKM theory, and they also reproduce the experimental selectivity. However, the master equation (which is reviewed in Section 5.1) was originally developed and used for studying the pressure and temperaturedependence of gasphase reactions, and when one tries to apply it in the liquid phase, there are huge uncertainties in what to use for the average energy transfer parameter and the collision frequency. Glowacki et al. varied these parameters to reproduce the experimental selectivity, and their adjusted collisional frequencies are around 3–10 ps^{−1}, and the adjusted 〈ΔE〉_{down} are 900–1100 cm^{−1}. We may compare this to the gas phase where the usually used value of 〈ΔE〉_{down} for small to medium sized molecules is around 100–500 cm^{−1}. These authors consider their adjusted energy transfer parameters “are broadly consistent with gas phase studies”, and their adjusted collisional frequencies are consistent with “experimental studies of solution phase intermolecular vibrational relaxation show fast intermolecular energy transfer to the solvent.”
Another very important question is the accuracy of pressuredependent VTST calculations with SSQRRK theory. Rigorously speaking, the pressuredependence calculation is beyond transition state theory; transition state theory just provides the key ingredients (highpressurelimit rate constant, and the effectively transformed energyresolved microcanonical SSQRRK rate constants) that are needed in the computation of pressuredependent rate constant.
One important assumption about the nonstatistical effects in such calculations is the model we adopted for treating the energy transfer. The collision efficiency we used (i.e., eqn (196) in Section 5) is derived from the solution of master equation with an exponentialdown model for vibrational energy transfer; the probability of an energy transfer collision reducing the energy from E to E′ is assumed to have an exponential functional form, with the larger the energy being transferred, the smaller the transfer probability. This has often been used satisfactorily, especially if one takes the uncertainty of the energy transfer parameter itself into consideration; a test for such a model on the simple CHF_{3} dissociation reaction^{198} (the dissociating products are HF and CF_{2}) and the C_{2}F_{4} dissociation reaction (which is a uphill from reactant to two CF_{2} on a PES without saddle point) at combustion temperatures and over a wide pressure range was consistent with the expected accuracy. Another deficiency in the current SSQRRK treatment is that we did not explicitly consider the J (angular momentum) dependence of the microcanonical rate constant, since the QRRK model is based on vibrational energy randomization; but the necessity of introducing the Jdependence in such a SSQRRK treatment needs further investigation since: (a) the appropriately adjusted energy transfer parameter (which is essentially an empirical parameter in most common practices) is able to compensate such dependence in the final falloff behavior; (b) SSQRRK is an efficient model to approximately capture all the physical effects contained in highpressurelimit thermal rate constants and transfer them into microcanonical rate constants, and the resulting microcanonical rate constants therefore implicitly also contain these physical effects and information about J.
7. VTST in condensed phases
7.1 Kinetics in liquid solutions
The central concept for reactions in liquid solutions is the free energy surface (FES), which is another name for the potential of mean force (PMF).^{180,263,264} The degrees of freedom for species that are of explicit interest (such as reactants) are called the primary degrees of freedom (primary system), and the remaining degrees of freedom (such as solvents) are called secondary degrees of freedom (secondary subsystem or environment). The PMF is a free energy function for the secondary subsystem, but an appropriately ensembleaveraged potential energy surface for the primary system.
In practice, the solvent degrees of freedom can be treated using an explicit solvation model, an implicit solvation model, or an explicit–implicit hybrid treatment.^{265,266} In the case of explicit solvation, the solvent molecules are treated at the full atomistic level; the configuration space has to be statistically sampled (via Monte Carlo or molecular dynamics) to provide an ensemble average. In an implicit solvation model, the solute molecules are embedded in the reaction field, which is the electrostatic field produced by the solvent when it is polarized by the solute; there are many reviews that cover implicit solvation models.^{267–272} In the explicit–implicit hybrid treatment or the socalled mixed discretecontinuum models, one explicitly adds a few solvent molecules (first solvation shell) around the solute molecule and then embeds them in a reaction field due to the remaining environmental molecules;^{273} it has been found that,^{274–276} in the case of strong interaction between solvent and solute (such as hydrogen bond, and molecular anions), adding a few explicit solvent molecules is able to yield better results as compared to using purely implicit solvation treatment.
Let R denote the 3N − 6 primary system coordinates, and P denote their conjugate momenta; the remaining degrees of freedom of the whole system are denoted by r for internal coordinates and p for conjugate momenta. The free energy G(S) of species S is

 (214) 
where the volume integral is carried out over the range of
R corresponding to species S; and
C is a geometryindependent constant. We can carry out this integral in two steps: first

 (215) 
and then

 (216) 
where
δ is the Dirac delta function, and
W(
R) is the PMF. The reason that
W(
R) is called the potential of mean force is that, when we take the derivative of
W(
R) with respect to
R, we have

 (217) 
where
U(
R,
r) is the total potential energy of the entire system, and the kinetic energy has been separated from the total Hamiltonian and has been integrated out over momenta. This result indicates that −∂
W(
R)/∂
R is the mean force acting on coordinate
R that is averaged over all other coordinates.
In the implicit solvation models, the PMF (the free energy surface) W(R) is

W(R) = V(R) + ΔG_{S}*(R)  (218) 
where
R denotes the 3
N − 6 internal solute coordinate,
V(
R) is the gasphase Born–Oppenheimer electronic potential energy surface, and Δ
G_{S}*(
R) is the fixedconcentration Gibbs free energy of solvation; it can be interpreted as the work of coupling a solute that is fixed in position to the equilibrated solvent at fixed temperature and pressure in an infinitely dilute solution.
^{277} The superscript “*” denotes a Gibbs free energy of solvation that corresponds to the same concentration (
e.g., 1 mol L
^{−1}) of the molecule both in the gas phase and in the liquidphase solution; that is an ideal gas at a gasphase concentration of 1 mol L
^{−1}, dissolving as an ideal solution at a liquid phase concentration of 1 mol L
^{−1}. In a practical calculations that employ an implicit solvation model (for instance SMD
^{278}), Δ
G_{S}* has already been included in the SCF singlepoint energy,
i.e., what we have calculated directly from a selfconsistent reaction field (SCRF) calculation is the free energy surface.
The molar Gibbs free energy of a species i in liquidphase solution is:^{277}

G_{i} = G_{i}* − TS_{lib,i}  (219) 
where
G_{i}* is the chemical potential of a species whose center of mass is constrained to a fixed position in the solution (it is called as the pseudochemical potential, and it is associated with the internal free energy of the species in solution including its coupling to the solvent); one may approximate
G_{i}* in as:

G_{i}* = U_{0} + G_{int}  (220) 
where
G_{int} is the internal thermal Gibbs free energy summed over all the configurations;
U_{0} is the zeropoint inclusive free energy surface, given in the harmonic approximation by:

 (221) 
where
F_{vib} is the number of vibrational degrees of freedom, and
ω_{m} is the vibrational frequency of mode
m. The entropy term
S_{lib,i} in
eqn (219) is the entropy of liberation of species
i, which is given by:

 (222) 
where
R is the gas constant,
c_{i} is concentration of
i in moles per unit volume,
N_{A} is Avogadro's constant. The last term of
eqn (219) is the free energy of liberation, which is formally identical to the translational molar free energy in an ideal gas.
In order to obtain the thermodynamic standardstate Gibbs free energy (which is directly related to partition functions and standard equilibrium constants), the standardstate Gibbs free energy of species B in solution (i.e., the chemical potential of B) should be computed as:

 (223) 
where
is the standardstate gasphase Gibbs free energy of species B, which is usually computed based on gasphase computed frequencies of the gasphase optimized geometries; but one can also use solutionphase computed frequencies of solutionphase optimized geometries.
^{279} The superscript “°” in
eqn (223) denotes the usual thermodynamic standard state, which is defined as follows: (a) gasphase molecules: all the molecules are ideal gases at a pressure of 1 bar, which corresponds to 0.0404 mol L
^{−1} at 298 K; (b) solutes in liquidphase solution: 1 mol L
^{−1}; (c) solvents: calculate from the density,
e.g., 55.5 mol L
^{−1} for water at 298 K. The quantity
in
eqn (223) is the amount that one must add to the fixedconcentration Gibbs free energy of solvation to change the standard state for the left hand side of the solvation process (the gas phase) to a 1 bar standard state from the 1 mol L
^{−1} standard state that would correspond to fixed concentration, while keeping the standard state for the righthand side of the solvation process (the liquidphase solution) at 1 mol L
^{−1}; the value of
is +1.9 kcal mol
^{−1} for solutes at 298 K, and it is +2.4 kcal mol
^{−1} for water (as a pure solvent) at 298 K.
The inclusion of tunneling in condensedphase calculations with explicit solvent does not require additional dynamical approximations although sampling raises new issues, as discussed in Section 7.3. With implicit solvent though, one only has the PMF, which is a statistical quantity, whereas in principle one should carry out tunneling on the true potential energy function. In other words, one should calculate the tunneling using the potential energy surface and average the dynamics, but with the PMF one has averaged the potential, and it is an additional approximation to carry out tunneling on the averaged potential. This additional approximation is called the canonical mean shape (CMS) approximation.^{280} With this approximation, VTST calculations are readily interfaced with various solvation model calculations.^{21,23,273,281,282}
Further complications could arise due to strong interactions between solute and solvent, and the dynamic solvent effect becomes nonnegligible. The key issue here is participation of the solvent in the reaction coordinate.^{263,283} With explicit solvent the reaction coordinate can have a component of solvent motion, but this is not possible when one uses the PMF. We should keep in mind that specifying the transition state dividing surface is the same as specifying the reaction coordinate since the reaction coordinate is the degree of freedom normal to this surface (the degree of freedom missing in this surface). If the best dividing surface has a reaction coordinate that includes a component of solvent motion, then a reaction coordinate without solvent motion must correspond to a dividing surface that leads to a larger transition state theory rate. Therefore the effect of including solvent motion in the reaction coordinate can be mimicked by adding a friction term that slows down the rate. Such friction terms are present in Kramers theory,^{284} Grote–Hynes theory,^{26,285–287} and other stochastic dynamics theories.^{288}
There are three approaches for computing the reaction rate constant for a liquidphase reaction:^{281,289} (1) separable equilibrium solvation (SES); (2) equilibrium solvation path (ESP); (3) nonequilibrium solvation (NES). Both SES and ESP are based on the equilibrium solvation assumption, which assumes that the solvation is in equilibrium with solvent at all points along the reaction path. In NES, on the other hand, the solvent coordinates are explicitly required for defining the reaction path, and this is called dynamic or nonadiabatic solvation.
In SES, one first finds the reaction path in gasphase, and then adds the free energy of solvation (at a fixed solute geometry) to the gasphase free energy for each point on the gasphase reaction path, using the gasphase optimized geometries and gasphase MEP. In ESP, the reaction path is determined by using the potential of mean force, which has been described in eqn (218); the position of the variational transition state in ESP not only moves along the gasphase reaction path but also perpendicular to it in the remaining 3N − 7 degrees of freedom; thus the ESP method is more flexible and more realistic than the SES method. Although ESP is a more complete treatment as compared to SES, the computationally less expensive SES often reasonably agrees with ESP^{281} (within a factor of 2, when the same solvation model^{269,273} is employed and quantum tunneling is considered). In NES, the specification of the reaction path explicitly involves one or more solvent coordinate as well as the solute coordinates.
Nonequilibrium solvation can be considered via a linearresponse theory for solvent friction effects, which can be incorporated with VTST and tunneling;^{290,291} the key concept in such linearresponse treatment, is the division of the system into an explicit subspace and an implicit bath, and a solvent coordinate (which is chosen based on physical basis and it is often preferable to use one or more collective solvent coordinates, that is a coordinate representing the solvent electric polarization field) is introduced in the Hamiltonian using the approximation of a single harmonic oscillator linearly coupled to the reaction coordinate.
7.2 Solid–vapor interface kinetics
Classical transition state theory has been applied to solids; we will limit our discussion to computing the migration and diffusion of atomic species on crystalline surfaces. Voter and Doll^{292} compared classical trajectory simulations results with classical TST results for diffusion constants. In trajectory simulations, the diffusion constants D are obtained by: 
 (224) 
where d is the dimensionality of the system; it is 2 for most surfaces and 1 for channeled surfaces like the (211) surface of bcc crystals, and 〈Δr^{2}(t)〉 is the meansquared displacement of the migrating atom, which can be obtained by averaging over trajectories. In practice, D is obtained by plotting 〈Δr^{2}(t)〉 with respect to time t, and fitting a line to the limiting slope. The temperature dependence of the diffusion constant is approximately described by an Arrheniustype equation: 
D = D_{0}e^{−Ea/kBT}  (225) 
If the motions of the migrating atoms are assumed to be independent random hops between adjacent binding sites, the diffusion constants can be computed by transition state theory as

 (226) 
where
l is the distance (hop length) between adjacent minimumenergy binding sites, and
k_{hop} is the hopping rate. Onedimensional harmonic conventional classical transition state theory yields

k_{hop} = n_{p}v_{0}e^{−V‡/kBT}  (227) 
where
V^{‡} is the energy difference between saddle point (between two binding sites) and local binding site minimum;
v_{0} is the harmonic frequency of the migrating atom moving along a pseudoonedimensional potential, in which all the other atoms are adiabatically relaxed; and
n_{p} is the number of binding sites accessible by a single hop. If we equate
V^{‡} to
E_{a} and equate
eqn (225) and (226), the parameter
D_{0} in the Arrheniustype equation for the diffusion constant is:

 (228) 
Transition state theory can also be used to compute the classical TST rate constant for transition from state A to state B without invoking the harmonic approximation:

 (229) 
where
f(
R) = 0 defines the dividing surface, and the ensemble average 〈
δ[
f(
R)]∇
f〉
_{A} is taken over the configuration space belonging to state A (TST).
^{293,294} Voter and Doll computed the ensemble average by the Monte Carlo method using the Metropolis algorithm as
^{295} 
〈δ[f(R)]∇f〉_{A} = f_{B}/ω  (230) 
where
ω is the width of the simulation box, and
f_{B} is the fraction of the Monte Carlo steps that lie inside the box.
VTST can be used, instead of conventional TST, to compute the hopping rate k_{hop}, and therefore it can give a better prediction for diffusion constants. The practical computations can be carried out with an embeddedcluster model.^{296–300} In such model systems, there is only one nonmetal atom. For instance, in the treatment of H atom diffusion on finite copper(100) planes,^{297,298} the PES is approximated as the sum of finite range H–Cu and Cu–Cu pair potentials. The H atom and the near six Cu atoms are allowed to move in the calculation, and other Cu atoms are fixed at experimental bulk geometry. The areas of these lattice planes and their number are taken large enough so that any atoms within the finite ranges of the interaction potentials of the seven nonfixed atoms (i.e., H and the six near Cu atoms) are included. The reactant is H located in one site (hollow formed by four Cu atoms) and the product is H in the nextdoor site; and the transition state structure is a bridge site between these two fourfold sites. Normal mode analysis is carried out along the minimum energy path, and canonical variational transition state is determined by maximizing the Gibbs free energy of activation; multidimensional quantum mechanical tunneling can also be included in the computation of the hopping rate to study the diffusion of hydrogen on metal surface.^{301–304}
For hydrogen diffusion on the Ni(100) surface, it has been experimentally observed^{305–307} that there exists a transition temperature, and for temperatures that are lower than the transition temperature, the diffusion coefficients were interpreted experimentally as nearly temperatureindependent due to tunneling. The quantitative definition of the transition temperature is the temperature at which the Arrhenius plot has maximum curvature. This zerotemperaturelimit, where the rate coefficient becomes a constant, has also been observed in organic reactions.^{308} However, the transition temperature was misinterpreted in the physics literature. Reanalysis^{300} of the lowtemperature situation shows that the experimental “transition temperature” actually signifies the transition from groundstatedominated tunneling to thermally activated tunneling. Thermally activated tunneling continues to dominate the rate constant at temperatures far above the transition temperature. In lowtemperaturelimit, only the vibrational ground state is populated, and the rate constant is constant because it simply corresponds to the rate of reaction by tunneling out of a single state. Under the harmonicoscillator approximation, the lowtemperaturelimit diffusion coefficients are:

 (231) 
where
l is the lateral distance between the two minimumenergy sites,
d is the dimensionality of the systems (which is usually 2 for surface diffusion),
σ is a symmetry factor that accounts for number of equivalent paths from reactants to products for the reaction (which is 4 for the (100) surface),
c is the speed of light,
^{R} is the vibrational frequency in wavenumber for the lowestfrequency hydrogenic vibration at the reactant well (which corresponds to the reaction coordinate over most of the reaction path away from the reactant well), and
P^{G}(
E^{R}_{0}) is the quantum transmission probability at the harmonic vibrational groundstate energy of the reactant well. At higher, but still low temperatures, where other states start to contribute, one can approximate the tunneling as occurring at a sequence of discrete energies in the reactant well, rather than out of a translational continuum.
^{296}
7.3 Enzyme kinetics
Unlike in the solution phase, where the solute and solvent can usually be clearly separated, in enzymes the partition of the entire system into primary and secondary is somewhat arbitrary. In practical computations of enzyme kinetics, one often uses a combination of molecular mechanics (MM) and quantum mechanics (QM), i.e., the QM/MM method. Enzyme–solvent systems can have a huge number of possible configurations, and these configurations make significant contributions to the free energy of activation; the configurations are usually sampled by MC or MD techniques. Ensembleaveraged VTST^{309} (EAVTST) with multidimensional tunneling was developed for modeling kinetics in enzymes, and we will briefly review EAVTST here; the readers who are interested in the detailed theory are directed to other previous publications with a higher concentration of focus on this topic.^{30,33,36,118,310–314}
EAVTST calculations can be divided into three stages, and each of the stages corresponds to a different level of completeness of the dynamical treatment.
Stage 1.
(Step 1) In this step we perform classical molecular dynamics (via, for instance, umbrella sampling) along a predefined reaction coordinate z (which can also be called a distinguished reaction coordinate, with a good example being the difference between internuclear distances of the forming and breaking bonds, and a more sophisticated example being a collective coordinate that is defined by the energy gap between the valence bond states corresponding to the reactant and product states) to produce the classical potential of mean force W_{C}(z) (free energy profile, or – more technically – the generalized free energy of activation profile). Once the classical potential of mean force is available, the generalizedtransitionstate (GT) classical potential of mean force activation energy ΔW_{GT,C}(T,z = z*) at temperature T is computed by: 
 (232) 
which is the difference between the maximum of classical potential of mean force and the reactant classical potential of mean force. Here “C” means classical.
(Step 2) In the first step, ΔW_{GT,C}(T,z = z*) has already included classical free energy contributions associated with all degrees of freedom that are orthogonal to the reaction coordinate; at a generalized transition state, the reaction coordinate itself does not contribute to the free energy (while all the vibrational modes that are orthogonal to reaction coordinate do contribute), but at the reactant it does. Thus, in the 2nd step, we add the missing contribution for the vibrational reactioncoordinatemode free energies of the reactant to the computed classical ΔW_{GT,C}(T,z = z*); and the soobtained quantity is the generalized transition state classical free energy of activation given by

ΔG_{GT,C}(T,z = z*) = ΔW_{GT,C}(T,z = z*) − G_{C,F}(T,z = z_{R})  (233) 
where G_{C,F}(T,z = z_{R}) is the classical free energy contribution of the reaction coordinate at its value at the reactant state z_{R}; it can be computed by calculating the free energy difference with and without this coordinate by not projecting and projecting the reaction coordinate from the Hessian matrix respectively.
(Step 3) The last part of this stage is to replace classical vibrational partition functions by quantum mechanical vibrational partition function, in order to include the quantization effects in vibrational free energies. In practice, this quantization is only done for an N_{1}atom subsystem; in particular, only 3N_{1} − 7 modes are quantized at the transition state, and 3N_{1} − 6 modes are quantized at the reactant. The Gibbs free energy of activation obtained this way is called the singlereactioncoordinate (SRC) quasiclassical (QC) generalizedtransitionstate (GT) free energy of activation, which is denoted as ΔG^{SRC}_{GT,QC}. The VTST rate constant for stage 1, which is denoted as k^{(1)}, for a unimolecular enzymatic reaction (EP → ES, or EP → E + P, where E, S, and P denote, respectively, the enzyme, substrate, and product, and ES is a Michaelis adduct) is thus:

 (234) 
Stage 2.
In this stage, we calculate ensembleaveraged recrossing and tunneling transmission coefficients in order to obtain the ensembleaveraged quasiclassical variational transition state theory rate constant. For EAVTST with tunneling (abbreviated as “T”), the rate constant is: 
k^{EAVTST/T} = γk^{(1)}  (235) 
where k^{(1)} is the VTST rate constant obtained in stage 1, and γ is ensembleaveraged transmission coefficient 
 (236) 
where M is the total number of ensemble members included in the average, Γ_{i} is the recrossing transmission coefficient of member i, and κ^{T}_{i} is the tunneling transmission coefficient of member i. In principle, M is large enough to converge the summation; in practice one can use M of about 20. During stage 2, the system evolves in a fixed field of its surroundings (these surroundings are called the secondary zone, i.e., solvent and part of the substrate and coenzyme); this is called the static secondaryzone approximation. This stage involves semiclassical reactionpath calculations for a subsystem (which is called the primary zone) within the secondaryzone being in frozen configurations that are taken from a quasiclassical transition state ensemble in stage 1. Notice that, in eqn (236), averaging the transmission coefficient over M ensemble members is equivalent to letting the secondary zone participate in the reaction coordinate.
Stage 3.
Usually, the static secondaryzone (also called “the frozen bath”) assumption works very well, and the nonequilibrium solvation effects are negligible; it is sufficient to stop at stage 2. As a further improvement, the entropic contributions from the secondaryzone can be included in the transmission coefficient. Then in the equilibrium secondaryzone approximation (which can be viewed as an analogous to a liquidphase kinetics calculation with the equilibrium solvation model), the free energy (potential of mean force) is further corrected by the bath free energy (which is obtained by equilibrating the secondaryzone along the reaction path of stage 2, and the bath free energy as a function of reaction coordinate is obtained via free energy perturbation simulations), and the recrossing transmission coefficient is computed based on this corrected free energy surface.
8. Selected applications
In applications, one requires a potential energy surface. There are two options. One option is to employ an analytic function, which may be obtained by a simple theory, such as the London equation,^{315} or by a fit to electronic structure calculations.^{316–324} The second option is direct dynamics.^{325,326} In direct dynamics “all required energies and forces for each geometry that is important for evaluating dynamical properties are obtained directly from electronic structure calculations.”^{71} In modern work, one almost always uses direct dynamics.
8.1 Combustion chemistry
VTST is a widely used tool for predicting the rate constants of elementary reactions that are important in fuel combustion, which is invaluable for selecting potential biofuels and fuel additives, for optimizing the combustion process, and for designing nextgeneration combustion engines. It has been demonstrated that anharmonicity^{150,327–330} plays an important role in controlling the VTST rate constants of hydrogen abstraction reactions. One interesting study that uses MS and MPVTST in analyzing the effect of hydrogenbonded transition states on rate constants (for gasphase hydrogen abstraction reactions) shows that,^{150} the conventional thinking of Hbonded TSs increasing rate constants, which is solely based on the energetic effects of the hydrogen bond that excludes the entropic effects, is not reliable; a hydrogen bond can reduce the entropy and thereby increase the free energy of the transition state (this is enthalpy–entropy compensation), and hence, a strong hydrogen bond interaction may lead to a slower reaction rate, and reliable conclusions about rate constants must be based on free energies of activation, not barrier heights or enthalpies of activation.
Among the many reactions that are important for ignition are the hydrogen abstraction reactions from the fuel molecules by various small radicals, such as H,^{78,88,331–333} HO_{2},^{334–337} OH,^{338–355} and O radicals.^{356–362} VTST can also supply rate constants for newly proposed reaction mechanisms for oxidation chemistry of fuels.^{363–370} VTST/MT has also been used to study the kinetics of isomerization reactions between organic radicals^{146,371–377} that are abundant in combustion flames. Guan and coworkers have also applied MSVTST to study hydrogen transfer between dimethyl ether and the methoxy radical,^{378} and combustion modeling of dimethyl ether.^{379} Kinetics of the unimolecular reactions of ethoxyethylperoxy radicals, which are main intermediates in the oxidation of diethyl ether under enginerelevant conditions, have also been studied using VTST.^{380} Klippenstein and coworkers have combined VTST with the master equation approach to study some of the key decomposition and association reactions in combustion.^{381–384}
8.2 Atmospheric chemistry
Studying atmospherically important reactions is critical to our understanding of climate change and air pollution. VTST together with pressuredependence theory is indispensible in providing accurate rate constants for global climate modeling. We have studied the chemistry of Criegee intermediates,^{385} whose importance in atmospheric science is widely appreciated, and understanding their fate is a prerequisite to modeling climatecontrolling atmospheric aerosol formation.
We have modeled the dissociation of C_{2}F_{4} with VRCVTST and SSQRRK theory;^{161} this represents a class of reactions and serves as an example for studying hydrofluorocarbon or Freon decomposition. The decomposition of hydrofluorocarbons or Freon is related to the destruction of the ozone layer.
We have also studied the reaction of SO_{2} with OH in the atmosphere;^{232} this reaction plays an essential role in acid rain and it is still attractive to experimentalists uptodate,^{386} and our temperature and pressuredependent rate constant computations help to resolve the discrepancies between various experimental studies. Cartoni et al. have studied the formation of HSO_{2}^{+} ion in the gas phase, which is of special interest because of an interesting reverse temperaturedependent kinetics.^{387} In the study of the important roles played by water and ammonia in promoting atmospheric reaction in atmosphere: the rate constant of watercatalyzed H_{2}SO_{4} + OH reaction is about 1000 times larger than that of the same reaction without water as reported by Long et al.;^{388} it has also been found that ammonia can accelerate this reaction by enhancing quantum tunneling.^{389} Using VTST with multidimensional tunneling, the unusual temperature dependence of the atmospherically important reaction OH + H_{2}S reaction has been understood, and predictions have been made for this reaction at higher temperatures.^{390} Loerting and coworkers have studied the proton transfer in small cyclic water clusters, which are important in understanding watercluster involved atmospheric reactions.^{391}
VTST is a powerful tool for filling the gap between limited experimental measurements and the need for rate constants in climate modeling at different altitudes, and in understanding previously proposed important atmospheric reaction mechanisms,^{392–399} and to provide rate constants that cover the tropospheric and stratospheric temperature and pressure range. As an example for using VTST/MT in modeling complex reaction mechanisms in atmospheric science and geoscience, Molina and coworkers^{400} have studied the heterogeneous reactions for chlorine nitrate hydrolysis on waterice surface in polar stratospheric cloud; their study incorporates proposed reaction mechanisms with reliable kinetics calculations, which lead to a good agreement of the estimated reaction probability with the experientially reported one for such a complex system.
8.3 Plasma chemistry
Siliconbased chemically active nanodusty plasmas are an important research subject in plasma physics and engineering.^{401–406} Engineers need thermochemical properties and rate constants to understand experimentally observed particle distributions, particle growth rates, and various transport phenomena in plasma. We have used density functional theory and MST to study the thermodynamic properties for various branched silicon clusters,^{407,408} and we found that the contribution from multiple conformational structures and torsional anharmonicity is significant; the conventional groupadditivity based estimations of thermodynamic properties were found to be unreliable. We have studied detailed chemical mechanisms for the growth of silicon hydride clusters,^{199,409} and we have computed the rate constants for all elementary steps by using VTST and multidimensional tunneling. We have also carried out a systematic benchmark study for testing various density functionals for their performance on elementary reactions in the polymerization of neutral silane with the silylene or silyl anion;^{199} this could serve as a reference for future siliconbased nanodusty plasma chemistry study. Oueslati et al. have compared the VTST/MT computed rate constants for hydrogen abstraction from tetramethylsilane,^{410} which can be viewed as a prototypical species in silanebased plasmas, and good agreements with experimental values was achieved.
8.4 Organic chemistry in solution
VTST with stateoftheart solvation models can be used to calculate reaction rates and elucidate reaction mechanisms in liquidphase solutions. As a recent example, we have used SMD implicit solvation model to study the E1cb mechanism for elimination reaction in liquid phase.^{264} The paper reporting these results contains an explanation of free energy surfaces, and using these we are able to distinguish between the concerted secondorder mechanism for β eliminations and nonconcerted mechanisms with discrete carbanion intermediates; distinguishing these experimentally is very difficult.
Liquidphase VTST calculations have also been carried out to investigate hydrogen migration in carbenes,^{411} hydrogen addition on heterocyclic ring,^{412} and hydrogen donation in reductive decarboxylation reactions.^{413} In the study of the kinetics of hydrogen transfer reaction for phenolic compounds in water,^{414} we have found that the M05 exchange–correlation functional^{415} yields more accurate rate constants than does M08HX, which is a density functional that is generally more accurate for gasphase reactions; this is possibly because the error introduced in the solvation model partially cancels the error in the lessaccurate density functional.
When using VTST to compute rate constants in the liquid phase, the accuracy of the final computed rate constants depends on many factors, and the accuracy of the solvation model is one of the most critical factors.
8.5 Enzyme reactions
Ensembleaveraged VTST with tunneling has been applied to study a number of enzyme reactions, including reactions catalyzed by yeast enolase,^{416} triosephosphate isomerase,^{417} methylamine dehydrogenase,^{418,419} alcohol dehydrogenase,^{309} thermophilic alcohol dehydrogenase,^{420} haloalkane dehalogenase,^{421} dihydrofolate reductase from E. coli,^{422,423} hyperthermophilic DHFR from Thermotoga maritima,^{424} xylose isomerase,^{425,426} shortchain acylCoA dehydrogenase,^{427} catechol omethyltransferase,^{428} glyoxalase I,^{429} coenzyme B12,^{430,431} HIV1 protease,^{432} 4oxalocrotonate tautomerase enzyme,^{433} hydride transfer between Anabaena Tyr303Ser FNR_{rd}/FNR_{ox} and NADP^{+}/H,^{434} hydride transfer from NADH to FMN in morphinone reductase.^{435} Some of this work is reviewed elsewhere.^{30,33,36,310,311,314,436}
One of the most important quantities in enzyme kinetics is the kinetic isotope effect (KIE), which is often used to determine mechanisms (primarily by helping to identify the slow step) and which often provides evidence for quantum mechanical tunneling. The ability for VTST to accurately computing such KIEs is strong evidence for the soundness of the theory.^{31,36,419,425,437–445}
One interesting finding is that, it is possible for the isotope D to tunnel more than H,^{446,447} and this can only be explained by multidimensional tunneling^{239,448} (as opposed to onedimensional tunneling) because if both isotopes tunnel along the same path with the same effective potential, H will always tunnel more. However, in multidimensional tunneling, both the tunneling paths and the effective potentials depend on all the masses in the system.
9. Summary
In this review, we have discussed the fundamentals of variational transition state theory, including multidimensional tunneling and methods for incorporating multiplestructure and torsionalpotential anharmonicity in VTST. We reviewed variablereactioncoordinate VTST as well as reactionpath VTST. We have discussed the recent extension of VTST for the convenient calculation of pressuredependent rate constants and falloff effects. We also briefly reviewed the use of VTST in condensed phases, including solid–gas interfaces, liquid solutions, and enzyme kinetics. A number of recent applications of VTST with its various recent developments are briefly reviewed as well.
Here, we summarize the major steps in kinetics studies. Please notice that for the sake of brevity, the following guide does not include reaction systems that may require special theoretical treatment.
Kinetics in a few words 
1 
Benchmark electronic structure methods (usually various density functionals) against reaction energies and barrier heights for your system. A good summary of the accuracy and the computational cost of hundreds of commonly used electronic structure methods for kinetics is given in ref. 415. 
2 
For reactions with barriers, use reactionpath variational transition state theory (RPVTST) with multidimensional tunneling (MT); for barrierless reactions, use variable reaction coordinate variational transition state theory (VRCVTST). Carry out direct dynamics rate constant calculations using one of these methods. 