Variational transition state theory: theoretical framework and recent developments

Junwei Lucas Bao * and Donald G. Truhlar *
Department of Chemistry, Chemical Theory Center, and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, USA. E-mail: junweilucasbao@gmail.com; truhlar@umn.edu

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 (MS-VTST), multi-path VTST (MP-VTST), both reaction-path VTST (RP-VTST) and variable reaction coordinate VTST (VRC-VTST), system-specific quantum Rice–Ramsperger–Kassel theory (SS-QRRK) for predicting pressure-dependent rate constants, and VTST in the solid phase, liquid phase, and enzymes. We also provide some perspectives regarding the general applicability of VTST.


image file: c7cs00602k-p1.tif

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) multi-structural and multi-path variational transition state theory with multidimensional tunneling, theory of pressure-dependent rate constants, torsional anharmonicity, and their applications in combustion and atmospheric chemistry; (2) multi-configuration pair-density functional theory (MC-PDFT), including its separated-pair (SP) version with the correlating-participating-orbital (CPO) scheme, and their applications in modeling transition metal compounds and spin-splittings of divalent radicals.

image file: c7cs00602k-p2.tif

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 calculations7–9 has provided a solid quantum mechanical foundation for transition state theory in terms of short-lived 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 multiple-surface 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 state-specific reaction rates or for predicting product internal-state 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 theory14–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 gas-phase reactions; this section uses a reaction coordinate associated with a reaction path, and it includes tunneling. The next four sections continue the treatment of gas-phase 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 condensed-phase 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 gas-phase variational transition state theory

A detailed pedagogical review of the fundamental theory underlying gas-phase 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 gas-phase 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 reaction-coordinate 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. State-selected 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 local-equilibrium assumption in the reactant region is assumed valid for most gas-phase 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 high-pressure limit rate constants. In the high-pressure 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 pressure-dependent, “thermal” refers to a Boltzmann distribution of bath gases, but there is a pressure- and temperature-dependent non-Boltzmann energy distribution of reactants. Thus, in the fall-off region, rate constants are functions of both pressure and temperature. We will discuss pressure dependence in Section 5.

It is convenient to introduce mass-scaled coordinates,22,45 which may also be called isoinertial coordinates. With such a coordinate system, the motion of an N-atom system in three dimensions can be reduced to that of a single-particle with a reduced mass μ moving in 3N-dimensional 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 amu1/2 if one used the familiar mass-weighted coordinates46 of spectroscopists. The motions of all the particles are governed by a potential energy surface V({ri}) that depends on the set of N three-dimensional positions ri 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 3N-dimensional generalized mass-scaled 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 6N-dimensional phase space represents an N-atom system, and the motion of each point is governed by the following set of Hamilton equations:

 
image file: c7cs00602k-t1.tif(1)
where the classical Hamiltonian is
 
image file: c7cs00602k-t2.tif(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:

 
image file: c7cs00602k-t3.tif(3)
where is the 6N-dimensional divergence operator, v is a 6N-dimensional 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 NR, and it can be evaluated by integrating the density of phase points over the phase space of the reactant:

 
image file: c7cs00602k-t4.tif(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
 
image file: c7cs00602k-t5.tif(5)
Then we can use Gauss's theorem to rewrite the volume integral in eqn (5) as a surface integral,
 
image file: c7cs00602k-t6.tif(6)
where n is a unit vector normal to the surface S pointing out of the volume R, and S is a (6N − 1)-dimensional hypersurface that separates the reactant region from the product region; this hypersurface is called the dividing surface.

The local one-way flux of systems going from reactants to products through this dividing surface is

 
image file: c7cs00602k-t7.tif(7)
where we consider the portion, S+, of the dividing surface for which the flux is positive. We define the coordinate q3N 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 one-way flux is evaluated as:
 
image file: c7cs00602k-t8.tif(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

 
ρ = ρ0eH/kBT(9)
where H is the Hamiltonian, kB is Boltzmann constant, T is temperature, and ρ0 is a constant. Plugging eqn (9) into eqn (8) and (4) yields
 
image file: c7cs00602k-t9.tif(10)
 
image file: c7cs00602k-t10.tif(11)

Now, we define a generalized transition state (GT) as a (6N − 2)-dimensional fixed-z hyperplane with the following Hamiltonian HGT:

 
image file: c7cs00602k-t11.tif(12)
The generalized-transition-state Hamiltonian HGT 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 HGT, i.e., HGT depends parametrically on z. Substituting eqn (12) into eqn (10), we obtain the one-way flux FGT(T,z*) through a generalized transition state fixed at z = z*, as
 
image file: c7cs00602k-t12.tif(13)
where GT indicates that the integral is carried out over the generalized transition state.

The local equilibrium one-way flux FGT(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 non-recrossing 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 FGT(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 so-called 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 kGTClassical(T), which is the flux per volume divided by the product of the concentration of reactants, obtained by the above procedure is

 
image file: c7cs00602k-t13.tif(14)
where V is the three-dimensional real-space 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
 
image file: c7cs00602k-t14.tif(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 = QClassical(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 VRP evaluated at z = z* on the reaction path (RP). The generalized transition state classical partition function per volume for the transition state becomes

 
image file: c7cs00602k-t15.tif(17)
The total Hamiltonian can be written as the sum of the Hamiltonians of the reactants, HA and HB, and using
 
NRX = ρX0h3NXXClassical(18)
where NRX is the number of atoms in X (NRA + NRB = N), and the constant ρ0 = ρA0ρB0, we obtain
 
image file: c7cs00602k-t16.tif(19)
in which X is A or B, X represents that the volume integral is carried out over reactant region, and HX is the Hamiltonian of the reactant X.

We then use eqn (15)–(18) to re-write eqn (14) as

 
image file: c7cs00602k-t17.tif(20)
Follow a similar derivation, the classical generalized TST rate constant for unimolecular reactions is:
 
image file: c7cs00602k-t18.tif(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 non-stationary point that is close to the saddle point on the reaction path is computed using:

 
x(s1 = ±δ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 normal-mode 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 non-stationary point x(sj) on MEP, then the next point can be computed by using the normalized gradient as:
 
image file: c7cs00602k-t19.tif(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 kGTClassical(T) since the local one-way flux provides an upper bound of the global reactive flux. The CVT rate constant is computed by

 
image file: c7cs00602k-t20.tif(24)

Therefore

 
image file: c7cs00602k-t21.tif(25)
which leads to the following CVT rate constant for bimolecular reactions
 
image file: c7cs00602k-t22.tif(26)
and for unimolecular reactions
 
image file: c7cs00602k-t23.tif(27)
Another way to understand the canonical variational approach is to re-write the classical generalized TST rate constant in terms of the standard-state Gibbs free energy of activation ΔGGT,0act(T,s).
 
image file: c7cs00602k-t24.tif(28)
where K0 is the reaction quotient (for the formation of TS) evaluated at the standard-state concentration of each species:
 
image file: c7cs00602k-t25.tif(29)
where, for gas-phase bimolecular reaction with unit cm3 molecule−1 s−1, (c°)−1 is kBT/p°, where p° is 1 bar; for liquid-phase reaction, the standard-state concentration is 1 mole L−1.

CVT is equivalent to maximizing the free energy of activation

 
image file: c7cs00602k-t26.tif(30)

In a microcanonical ensemble (NVE ensemble), the microcanonical energy-resolved classical rate constant (which is obtained by firstly integrating the E, J-resolved rate constant kGT(E,J,s) over the angular momentum J distribution) is

 
image file: c7cs00602k-t27.tif(31)
where NGT(E,s) is the phase-space 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:
 
image file: c7cs00602k-t28.tif(32)
and therefore the location of microcanonical variational transition state is energy-dependent. To obtain canonical μVT rate constant, one needs to perform the following integration:
 
image file: c7cs00602k-t29.tif(33)
where PR(E) is the normalized population distribution of the reactant:
 
image file: c7cs00602k-t30.tif(34)
where QR(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
 
kCVT(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 pz 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 3-dimensional 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 CH2), and a few polyatomics have large rotational constants, for which the classical approximation is not always sufficiently accurate. For H2, 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 minimum-energy structure, whereas real molecules often have multiple structures (conformations) due to torsions and/or rings. In some cases, the single-structure 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 harmonic-oscillator approximation truncates the series expansion of the potential energy at the second-order terms). Vibrational anharmonicity can be included by vibrational self-consistent-filed theory (VSCF), vibrational perturbation theory (VPT), vibrational configurational interaction (VCI) and other VSCF-based methods,57 but this can get expensive. To avoid the cost of those higher-level vibrational calculations, it is convenient use a scale factor58 (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 quasi-harmonic 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 non-separability in a separable-mode 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 factors58 are developed based on a database that consists of 15 representative small molecules, whose experimental harmonic vibrational frequencies, fundamental frequencies and zero-point 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 degeneracy-corrected,60 second-order61–63 vibrational perturbation theory (HDCVPT2); and it is equal to the ratio of HDCVPT2 computed anharmonic ZPE to the harmonic-oscillator 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 zero-point energy are primarily dictated by the stretching frequencies since they are larger than the bending and torsion frequencies and hence their contributions dominate the zero-point energy. Since stretches usually have negative anharmonicity, frequency-scaling factors are usually smaller than unity, which means that partition functions are increased by using a frequency-scaling 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 high-frequency mode/non-torsional 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 above-mentioned HDCVPT2 method.

The quantized total partition function of species X can be then written as

 
Q(X) = Qtrans(X)Qrot(X)Qvib(X)Qelec(X)(36)
with the vibrational partition function in the QHO approximation given by
 
image file: c7cs00602k-t31.tif(37)
where for the reactant and product, the vibrational degree of freedom F = 3N − 5 for linear molecule, or 3N − 6 for nonlinear molecule, 3N − 6 for linear transition states, and 3N − 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 VMEP(s). We also define the zero of energy for the electronic partition function as VMEP(s), and hence, without considering spin–orbit coupling or low-lying 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 low-lying 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, QGT 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

 
Qelec = g0 + g1eE1/kBT +⋯(38)
where g0 and g1 are the degeneracy of ground state and first excited state, and E1 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 ground-state degeneracy. Notable exceptions are reactions involving halogen atoms and the OH radical.

Next we need to consider the vibrational zero-point 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 zero-point 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 ground-state potential energy VGa as follows:

 
VGa = VMEP(s) + εG(s)(39)
where εG(s) is the zero-point 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 lowest-energy 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 ei/2kBT in the numerator of eqn (37). In the rate constant expressions, the barrier height is then ZPE-exclusive. 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 ZPE-inclusive, i.e., it is based on VGa.)

The quantized CVT rate constant (which we denote as CVT/T, where “/T” means tunneling) can be then written as

 
image file: c7cs00602k-t32.tif(40)
 
image file: c7cs00602k-t33.tif(41)
where κT is the (temperature-dependent) tunneling transmission coefficient, which will be discussed in Section 2.3. image file: c7cs00602k-t34.tif 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, image file: c7cs00602k-t35.tif is simply the potential energy of the variational transition state; QGT and QA are the translational-motion-excluded quantized partition functions; ΦR is the product of the reactants’ partition functions per volume, which is equal to
 
image file: c7cs00602k-t36.tif(42)
in which mAmB/mGT is just the reduced mass of relative translational motion of A and B (because mGT = mA + mB). An alternative (equivalent) way for computing partition functions is that, in eqn (40) and (41), the partition functions QGT and QA 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 re-written as:

 
kCVT/T(T) =κTΓCVTkTST(43)
where kTST 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 work245 in the very early days.) If we ignore tunneling (i.e., we set κT = 1), then the CVT variational transmission coefficient is
 
ΓCVT = kCVT/kTST(44)

The value of ΓCVT characterizes the importance of variational effects (recrossing) in the studied system, and it is smaller than unity. Under the harmonic-oscillator approximation, CVT variational transmission coefficient is computed as

 
image file: c7cs00602k-t37.tif(45)
where QVT-HO is the vibrational partition function of variational transition state, Q‡-HO is the vibrational partition function of the saddle point; VVT 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 temperature-dependent, and so does the ΓCVT.

Similarly, for microcanonical variational transition state theory (μVT), we have

 
image file: c7cs00602k-t38.tif(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 number79

 
image file: c7cs00602k-t39.tif(47)
and the corresponding VTST rate constants are expressed as (in which the rotational partition function excludes the rotational symmetry number)
 
image file: c7cs00602k-t40.tif(48)
 
image file: c7cs00602k-t41.tif(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 CH4 + OH → CH3 + H2O, the number of equivalent ways for OH to abstract the H atom from CH4 is 4 (CH4 has four chemically equivalent H atoms), but the reaction symmetry number for forward reaction is 12 (the order of rotational subgroup for Td point group (CH4) is 12, and rotational symmetry number for OH is 1 and for TS (Cs point group) is 1); but in the case of CH4 + O → CH3 + OH, because the point group of the TS is C3v (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 CH4. The reaction symmetry number should also include the contributions from the non-superimposable mirror images.80 However, if the MS-T method137,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 MS-T partition function and we should not double count them.

We should remind ourselves that after this non-rigorous 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 reaction-coordinate 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 reaction-coordinate 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 heavier-atom 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:
 
image file: c7cs00602k-t42.tif(50)
where E0 is the higher of either the ground-vibrational-state energy of the reactants or the ground-vibrational-state energy of the products (this is the lowest energy at which reaction can occur), PT(E) is the quantum tunneling probability, and PC(E) is the classical transmission probability. We assume that the effective potential for tunneling is the vibrationally adiabatic ground-state potential energy curve VGa of eqn (39), whose maximum value is called VAG. Then, if the transition state were at the maximum of the vibrationally adiabatic ground-state potential energy curve, we would have
 
image file: c7cs00602k-t43.tif(51)
where VAG is the maximum of the vibrationally adiabatic ground-state potential energy VGa along the reaction coordinate. But in general, the CVT transition state is not at the maximum of the vibrationally adiabatic ground-state potential energy; thus we have to replace VAG by image file: c7cs00602k-t44.tif where image file: c7cs00602k-t45.tif is the location of the CVT transition state at temperature T.

The above definition of the PC(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

 
image file: c7cs00602k-t46.tif(52)
As for the quantum tunneling probability PT(E), it satisfies:
 
image file: c7cs00602k-t47.tif(53)
The region below barrier is the tunneling region, and the region above barrier is non-classical reflection region. The remaining problem is approximating the quantum tunneling probability PT(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: one-dimensional 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 one-dimensional tunneling models do not. The one-dimensional models reviewed in this article include the Eckart tunneling model,81,82 the Wigner tunneling model,83 the one-dimensional truncated parabolic tunneling approximation;84 and multidimensional tunneling models include zero-curvature tunneling (ZCT),49,50 the parabolic model with an effective frequency,85 small-curvature tunneling (SCT),86 large-curvature tunneling (LCT),71,87–90 microcanonically optimized multidimensional tunneling (μOMT),71 and least-action tunneling (LAT).91–93 All of the one-dimensional tunneling models and the ZCT tunneling model ignore the reaction-path curvature; this is not a good assumption for many cases. The effects of the reaction-path 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 approximation94 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; uniformization95,96 is the process by which one include higher-order 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 one-dimensional models.


2.3.3.1 One-dimensional 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
 
image file: c7cs00602k-t48.tif(54)
where
 
y = exp[α(ss0)](55)
 
B = [Vmax1/2 + (Vmax − ΔV)1/2]2(56)
in which ΔV is the value of the potential at s = +∞ relative to s = −∞, Vmax is the value of the potential at its maximum relative to s = −∞. And the constant s0 is defined to give the maximum value of V(s), (i.e., Vmax) at s = 0:
 
s0 = (2α)−1[thin space (1/6-em)]ln(1 − ΔV/Vmax)(57)
The intrinsic barrier height Vb described by the potential is:
 
Vb = Vmax − 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
 
image file: c7cs00602k-t49.tif(59)
where the second derivative Fs, in practice, can only be approximately computed by:
 
Fs = 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 imaginary-frequency 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 Fs 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 PT(E) can be analytically evaluated, which provides computational simplicity, but it may yield quantitatively unsatisfactory results. The PT(E) given by Eckart potential is:
 
image file: c7cs00602k-t50.tif(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 zero-curvature 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 VMEP, 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 VMEP, which is probably why it often overestimates the tunneling probability compared to ZCT. (The VMEP curve tends to be thinner for the nearly symmetric barriers that have the greatest tunneling.)

Although PT(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

 
image file: c7cs00602k-t51.tif(65)
The corresponding tunneling probability, which can be derived from eqn (61), is
 
image file: c7cs00602k-t52.tif(66)
where
 
ε = E/Vmax(67)
 
image file: c7cs00602k-t53.tif(68)
If we assume the barrier is very high or very broad such that γ ≫ 1, then we have
 
image file: c7cs00602k-t54.tif(69)
Plug the above tunneling probability into eqn (52), we can get
 
image file: c7cs00602k-t55.tif(70)
in which
 
x = γ(EVmax)/Vmax(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,22i.e.,
 
image file: c7cs00602k-t56.tif(72)
and finally an analytical expression can be arrived:
 
image file: c7cs00602k-t57.tif(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 VMEP – not to VGa.

Eqn (73) can be further simplified since t/sin(t) ≈ 1 + t2/6 for t ≪ 1, we obtain the simple Wigner tunneling approximation:83

 
image file: c7cs00602k-t58.tif(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 E0 = −∞, 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|ν|/(2kBT) = π(75)
Skodje and Truhlar84 have eliminated this issue by using the uniformedized WKB semi-classical tunneling approximation, and they have further generalized the result to be applicable for any unsymmetrically truncated parabolic-type barrier that is defined as:
 
image file: c7cs00602k-t59.tif(76)
and the intrinsic barrier Vb is
 
Vb = Vmax − max[0,ΔV](77)
The Skodje–Truhlar one-dimensional truncated parabolic tunneling transmission coefficient κT(T) is given by:
 
image file: c7cs00602k-t60.tif(78)
which is a singularity-free generalization of the result given by eqn (73). The above results (eqn (76)–(78)) are for one-dimensional tunneling; later, Skodje et al. also showed that this parabolic model could also be generalized so that it includes the reaction-path curvature with the small-curvature 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 paper85 (which is valid only if the reaction-path 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 one-dimensional models for quantitative studies; and they were reviewed in this sub-section 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 mass-scaled coordinate system is one manifestation of the coupling.41,99 A major consequence of reaction-path 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 zero-curvature 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 VGa(s), which contains the vibrational ground-state energies of all the vibrational modes that are orthogonal to the reaction coordinate s, and this zero-point energy depends on s. The contribution of vibrations transverse to the reaction path to VGa(s) is the reason why the SCT approximation is multidimensional. The ZCT tunneling probability is given by:

 
image file: c7cs00602k-t61.tif(79)
where E0 is the maximum between the VGa of the reactants and of the products; s* is the location of the variational transition state, where the VGa reaches its maximum. The magnitude of the imaginary action integral is
 
image file: c7cs00602k-t62.tif(80)
where s< and s> are the two turning points, at which VGa = 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 reaction-path 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 small-curvature 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 centrifugal-dominant small curvature approximation), which is used in modern VTST calculations, μeff is determined by reaction-path curvature as follows:86,98

 
image file: c7cs00602k-t63.tif(81)
in which ā(s) = |κ(s)[t with combining macron](s)|, and the reaction-path curvature |κ(s)| is image file: c7cs00602k-t64.tif, where BmF(s) is the curvature component of vector BF(s), which represents the coupling between the reaction coordinate s and a mode m that is perpendicular to it:
 
image file: c7cs00602k-t65.tif(82)
where ni(s) is the component i of the unit vector perpendicular to the generalized transition state dividing surface at s and LGTi,m(s) is the i-th component of the eigenvector for vibrational mode m perpendicular to n(s) at s.

The effective harmonic potential V for the vibrational mode u1 (which is constructed by a local rotation of the vibrational axes such that BF(s) lies along with it, and by construction, the curvature coupling in all other vibrational coordinates u2 to uF−1 are zero in this rotated coordinate system) is

 
image file: c7cs00602k-t66.tif(83)
where ωm(s) is the vibrational frequency of mode m at reaction coordinate s. The associated turning point [t with combining macron](s) for zero-point motion in this harmonic potential is then:
 
image file: c7cs00602k-t67.tif(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 turning-point derivative terms dtm(s)/ds are evaluated at s by a two-point central-difference method using the turning-point 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

 
image file: c7cs00602k-t68.tif(85)
where δs is the step size of the reaction coordinate.

The corner-cutting model that led to the effective reduced mass μeff(s) does not apply to energies above the maximum of VGa, it is therefore inappropriate to compute the transmission probability at the non-classical reflection region based on PSCT(E). The transmission probability in the non-classical reflection region is computed by

 
PSCT(E) = 1 − PSCT(2VGaE), VGa(s*) < E ≤ 2VGa(s*) − E0(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 ground-state vibrationally adiabatic potential energy curve). The large-curvature tunneling approximation71,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 vibrational-ground-state 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 reaction-coordinate turning point. The classical reaction-coordinate turning points are the points at which the ground-state 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 reaction-path 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 reaction-path curvature, the tunneling probability PLCT(E) could be smaller than PSCT(E), because SCT is a more accurate treatment in this limit region. At intermediate reaction-path 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 PLCT(E) and PSCT(E). When one uses μOMT, one often finds some contribution form large-curvature tunneling at low energy, even if the more important higher-energy tunneling is dominated by small-curvature 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 least-action 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 multi-dimensional 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 non-converged 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 two-step reaction:
 
image file: c7cs00602k-t69.tif(87)
The first step is the formation of an intermediate (I), with a forward rate constant of k1; the intermediate I further proceeds to yield the product P, with a rate constant of k2. The overall reaction rate for the formation of product is:
 
image file: c7cs00602k-t70.tif(88)
and if the intermediate I is in equilibrium with A and B, we would have
 
image file: c7cs00602k-t71.tif(89)
where K is the equilibrium constant for the first step, and K = k1/k−1. And thus the overall bimolecular rate constant is
 
koverall = k2K(90)
With (variational) transition state theory, the rate constant k2 is
 
image file: c7cs00602k-t72.tif(91)
where κT2 is the tunneling transmission coefficient computed by using the intermediate I as the reactant. The equilibrium constant K of the first step is
 
image file: c7cs00602k-t73.tif(92)
and thus, the overall bimolecular rate constant is equivalent to
 
image file: c7cs00602k-t74.tif(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 κT2 in eqn (93) is replaced by κT1, 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 κT1 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 κT2 can be greater than unity. This shows that the well-known dictum that one can ignore an intermediate that occurs before the rate-determining 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 two-step unimolecular reaction). This is a very general issue since most bimolecular reactions have barrier, and a ground-electronic-state 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 high-pressure limit, there may be inadequate collisional stabilization of the complex. In the low-pressure limit, there are no stabilizing collisions, and therefore there are no tunneling contributions at energies below the ground-state energy of the reactants121 (which would be the lowest possible energy that a complex could have). The low-pressure limit and the high-pressure limit are easily handled as limiting cases, but at intermediate pressures one must use the theory of pressure-dependent 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 reaction-coordinate non-separability, namely (i) the already-mentioned reaction-coordinate 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 non-separability 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 one-way 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 reaction-coordinate motion.

In general, then there are two contributions to the transmission coefficient, and we write

 
γ = κΓ(94)
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
 
kCVT/T = κkCVT(95)
where kCVT is the quasiclassical VTST rate constant. However, sometimes we use a different notation, namely
 
kCVT/T = κΓCVTkTST = γkTST(96)
where kTST is the quasiclassical conventional VTST rate constant, and
 
ΓCVT = kCVT/kTST(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 no-recrossing 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 one-way 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 = kexact/kTST(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 multi-dimensional quantum tunneling transmission coefficients rely on the computed vibrationally adiabatic ground-state energy along the reaction path; the generalized normal mode analysis is performed at evenly spaced points along the minimum energy path. For non-stationary 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 reaction-coordinate degree of freedom). One should keep in mind the fact that the projected Hessians in this generalized normal mode analysis are coordinate-system-dependent,75,76 and so are the bound-mode vibrational frequencies of non-stationary points.

The reason that the vibrational frequencies along the MEP are coordinate-dependent 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 non-stationary points. These two different choices of the reaction coordinate are related by76,122

 
image file: c7cs00602k-t75.tif(99)
in which F is the number of internal coordinates, and F − 1 excludes the reaction coordinate s (which is also called qF). The qi represent curvilinear (internal) coordinates, which measure distortions away from the MEP, and by definition they vanish on the MEP. The bij involve second-order partial derivatives of s′ with respect to qi, and they are in general nonzero. The double summation runs over a complete, non-redundant set q1, q2,…, qF−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
 
image file: c7cs00602k-t76.tif(100)
in which q ≡ {q1,q2,…, qF−1,s} and q′ ≡ {q1,q2,…, qF−1,s′}. Clearly, at stationary points, where the first-order 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 non-stationary points, where the (∂V/∂s) term does not vanish, the frequencies are coordinate-system-dependent.

And thus, the important consequence is that, the results of the VTST/MT calculations depend on the choice of the coordinates75,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 path75,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 program42 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 ground-state energy (VGa) 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 above-mentioned 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 energy-dependent 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 multi-dimensional tunneling. In Polyrate, one could choose either a set of 3N − 6 (where N is total number of atoms) non-redundant internal coordinates75 (via the “CURV2” option), or redundant internal coordinates124 (via the “CURV3” option), as curvilinear coordinates for performing the generalized normal mode analysis. A set of physically meaningful and well-chosen internal coordinates is usually able to give an imaginary-frequency-free 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 re-orient the generalized-transition-state dividing surface (via the “RODS” option) to eliminate the imaginary frequencies along the path.125 As for the computations of the multi-dimensional 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 non-converged 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 State-selected 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 vibrational-state-selected rate constants.126,127 Liu recently reviewed the experimental progress made in studying the mode-selected reaction dynamics.128 In this section, we briefly review using VTST to study state-selected 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 state-selected mode, only one term contributes. Another is based on non-adiabaticity. 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 non-stationary points, F includes all the bound-mode vibrational degrees of freedom that are orthogonal to MEP, which excludes the reaction coordinate),

 
image file: c7cs00602k-t77.tif(101)
where m labels a vibrational mode, and nm is the number of quanta in that mode, and nmaxm is the number of quanta used to converge the summation. If the vibrational energy levels are computed by harmonic oscillator (HO) approximation
 
image file: c7cs00602k-t78.tif(102)
where vm 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 mode-selected VTST, one modifies eqn (101) to
 
image file: c7cs00602k-t79.tif(103)
in which the modes from m = 1 to m = Fr are the selected modes (i.e., each of these modes is restricted to have a specific number of excitation quanta equal to nrm); and the remaining FFr modes are treated as usual with the summation over all the vibrational quantum numbers nm 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 {nrm}. In adiabatic mode-selected 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 bound-mode vibrational frequencies along the reaction path. This same restriction is also used to generate state-selected 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 state-selected VTST is computed based on the following state-selected vibrationally adiabatic potential:

 
image file: c7cs00602k-t80.tif(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 reaction-path curvature; large reaction-path curvature causes strong vibrational nonadiabaticity,129 and in such a case, one should use the second way to treat state-selected reactions, namely the diabatic one126 with the vibrationally diabatic potential given by

 
image file: c7cs00602k-t81.tif(105)
which is similar in form to the expression for the vibrationally adiabatic potential Vga, but instead of using adiabatic frequencies, it uses diabatic frequencies which change smoothly along the reaction coordinate. When the adiabatic frequencies vm(s) show an avoided crossing in the frequencies vmvs. s diagram (this is called the vibrational-mode 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 + H2 → H2O + 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 highest-vibrational-frequency modes (which are denoted as m = 4, 5) with the same symmetry have an avoided crossing at s0 (which is about −0.16 Å). The diabatic frequencies vdm(s) are then defined as:
 
image file: c7cs00602k-t82.tif(106)
where θ is the step function, and δ is Kronecker delta.

State-selected 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 treatment131 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 reaction-path 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 lowest-energy structure and on the (quasi-)harmonic oscillator approximation. A harmonic oscillator has only one minimum-energy structure, but real species often have multiple structures, e.g., gauche and trans 1,2-dichloroethane. 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: multiple-structure anharmonicity and torsional potential anharmonicity. The multiple-structure 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 (MS-T) method, which has two versions: MS-T with uncoupled torsional potential method (MS-T(U) method), and MS-T with coupled torsional potential method (MS-T(C) method). We recommend the second version of this method, which is based on a coupled torsional potential.138 Including the MS-T treatment of anharmonicity in VTST will be called multistructural VTST (MS-VTST) or multipath VTST (MP-VTST), 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 MS-T partition function), we will generate nt 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 nt. All of the distinguishable conformers that are generated by torsions are included in the computation of the MS-T partition function, including non-superimposable 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 (MS-LH).139 However, treating low-frequency 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 high-temperature limit:

 
image file: c7cs00602k-t83.tif(107)
and this is the vibrational partition function of a classical harmonic oscillator.

At very low temperatures where kBT, an internal rotation can be reasonably well treated with the harmonic-oscillator approximation at each minimum of the torsional potential; at very high temperature where kBT is much larger than the internal rotational barrier W, and the internal rotation becomes a free rotor. At intermediate temperature where < kBT < W, torsional degrees of freedom should be treated as hindered rotors.

Here we review the MS-T method based on a coupled torsional potential (i.e., MS-T(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 MS-T rotational–vibrational (which can be abbreviated as rovibrational) partition function is computed as the summation of the single-structure torsional rovibrational partition function (SS-T) of the unique structure j. In the equations used for the MS-T 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 non-superimposable mirror image that can be generated via torsion also counts as a unique structure. The MS-T partition function is:

 
image file: c7cs00602k-t84.tif(108)
where the SS-T rovibrational partition function for distinguishable structure j is defined as
 
image file: c7cs00602k-t85.tif(109)
where Qrotj is the rotational partition function of structure j, QHOj is the harmonic-oscillator (or quasi-harmonic-oscillator) vibrational partition function of structure j, Uj is the relative Born–Oppenheimer equilibrium energy of structure j with respect to the Born–Oppenheimer equilibrium energy of the lowest-energy conformer (which is denoted as the global minimum or GM), t is the total number of torsions, and fj,η 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 single-structure harmonic-oscillator rovibrational partition function (SS-HO) 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:

 
image file: c7cs00602k-t86.tif(110)
where σrot,j is the overall rotational symmetry number for structure j (each conformer j has its own symmetry number), and Irotj,1, Irotj,2, and Irotj,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

 
image file: c7cs00602k-t87.tif(111)
where qRC(C)j,η is the reference classical partition function of coupled torsional motion η for structure j using the coupled effective torsional barrier, and qCHO(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 QFCC (where all t torsions and the overall rotation are coupled):

 
image file: c7cs00602k-t88.tif(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 QFCC: (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
 
image file: c7cs00602k-t89.tif(113)
where, in the coupled MS-T calculation, Wj,τ 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 Mj,τ 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 nearly-separable uncoupled –CH3 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, Mj,τ is generally a conformer-dependent non-integer, 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

 
image file: c7cs00602k-t90.tif(114)
where I0 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 Mj,τ values for each torsion. We divide the t torsions into two types: nearly separable (NS, for instance CH3 groups) and strongly coupled (SC). We use the notation NS:SC = tNS:tSC to denote that tNS torsions are treated as nearly separable and tSC 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 tSC 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

 
image file: c7cs00602k-t91.tif(115)
where tSC 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 Mj,τSC by considering each torsion separately; then we replace all Mj,τSC for strongly coupled torsions of a given conformer j by a single MSCj equal to
 
image file: c7cs00602k-t92.tif(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,141eqn (116) is computed by:
 
image file: c7cs00602k-t93.tif(117)
In this sampling algorithm, a total number of Ntot 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 Nj is the number of points that are assigned to structure j.

Here we give some numerical examples of the MSCj 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 MSCj values, because they are treated as nearly separable with M = 3. For 1-pentyl radical (for which the conformational search is carried out by M06-2X/MG3S), the MSCj values for the SC torsions in its eight distinguishable conformers (which excludes 7 distinguishable mirror images) are, from lowest-energy conformer to highest-energy 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 tert-butanol by HO2 radical (the products are (CH3)2C(OH)CH2˙ radical and H2O2) possesses five SC torsions and has 46 distinguishable conformers (i.e., 23 pairs of non-superimposable mirror images, which are found by M08-HX/MG3S); the MSCj values range from 2.00 to 2.45, with a mean value of 2.17 and a standard deviation of 0.13. For (S)-sec-butanol (for which the conformational search is carried out by M08-HX/MG3S), the MSCj values for the SC torsions in its nine distinguishable conformers (from the lowest-energy conformer to the highest-energy 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 qCHO(C)j,η with coupled frequencies is:

 
image file: c7cs00602k-t94.tif(118)
where F is the number of vibrational degrees of freedom, ωj,m is the vibrational frequency of normal mode m, and [small omega, Greek, tilde]j,m are the non-torsional vibrational frequencies; these frequencies are computed by diagonalizing respectively the full Wilson–Decius–Cross GF matrix and its non-torsional submatrix with non-redundant internal coordinates.

The final result is

 
image file: c7cs00602k-t95.tif(119)

In general, for a system with multiple coupled torsional degrees of freedom, the accurate partition function is unaffordable to compute. Nevertheless, for two-dimensional systems, i.e., molecules with two coupled torsional degrees of freedom, where the highly accurate partition functions are available via the extended two-dimensional torsion method144 (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 MS-T method with coupled torsional potential, i.e., MS-T(C) method, has been confirmed; and it is concluded that “If the E2DT method is not affordable, the MS-T(C) method is the best option to calculate rovibrational partition functions including torsional anharmonicity”.

To carry out MS-T 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 dual-level MS-T,145 which requires the user to carry out a conformational search with a chosen low-level theory (for instance, a semiempirical molecular orbital method) and then re-optimize the 30% energetically lowest low-level structures at a higher level. Dual-level MS-T is able to combine the information from a low-level theory with that from a higher-level theory to obtain an approximate MS-T partition function that agrees with the full high-level MS-T partition function within a factor 2 to 3; an even more encouraging feature is that we have shown that when one is using dual-level MS-T 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 dual-level method canceled out to a large extent), very satisfactory results (as compared to using full high-level MS-T partition functions) can be obtained.

3.2 Multistructural variational transition state theory

In multistructural variational transition state theory146,150 (MS-VTST), we replace the single structure harmonic oscillator partition functions in the expression of the VTST rate constant, by MS-T partition functions. And therefore, the multistructural and torsional anharmonicity is incorporated in the VTST theory.

We now derive the MS-VTST 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

 
image file: c7cs00602k-t96.tif(120)
where QGT-SSHO and QR-SSHO are respectively the single-structural (lowest-energy structure) harmonic-oscillator 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 re-write this expression in terms of the recrossing transmission coefficient and the conventional (single-structural) TST rate constant
 
image file: c7cs00602k-t97.tif(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
 
image file: c7cs00602k-t98.tif(122)
where Q‡-SSHO is the single-structural harmonic-oscillator rovibrational partition function for the lowest-energy saddle point. In order to take the MS-T effects into VTST calculation, we replace the single-structural transition state theory rate constant with the multistructural TST (MS-TST) rate constant in eqn (121). The multistructural TST rate constant is defined by replacing the single-structural partition functions with MS-T partition functions
 
image file: c7cs00602k-t99.tif(123)
which can be written in terms of the multistructural and torsional anharmonicity factor FMS-Tfwd of forward reaction:
 
kMS-TST(T) = FMS-TfwdkSS-TST(T)(124)
 
image file: c7cs00602k-t100.tif(125)

The multistructural canonical variational transition state theory (MS-CVT) rate constant is defined as the following:

 
kMS-CVT/T(T) = κTΓCVTkMS-TST(T) = FMS-TfwdkSS-CVT/T(T)(126)
where the single-structural CVT rate constant is
 
kSS-CVT/T(T) = κTΓCVTkSS-TST(T)(127)
Notice that, in the MS-CVT theory presented here, we approximate the ratio of the MS-T partition function to the SS-HO partition function at the variational transition state as the ratio at the conventional transition state; and the variational transmission coefficients are computed based on single-structural (quasi-) harmonic-oscillator 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 multiple-structural anharmonicity (which requires optimizations and generalized normal mode analysis for all the possible generalized conformers of the non-stationary points) along the reaction path, which is simply unaffordable for large or even medium-sized transition state structures. For bimolecular reactions, the factor FMS-Tfwd for forward reaction is defined as:
 
image file: c7cs00602k-t101.tif(128)

For the reverse reaction, we can also define the FMS-Trev 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 MS-T standard-state reaction equilibrium constant K°,MS-T being the product of FMS-Trxn factor for the reaction (which is equal to FMS-Tfwd/FMS-Trev), with the single-structural harmonic-oscillator standard-state reaction equilibrium constant K°,SS-HO (which can be derived from SS-HO standard Gibbs free energy of reaction image file: c7cs00602k-t102.tif); and the MS-T standard Gibbs free energy of reaction image file: c7cs00602k-t103.tif.

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. MS-VTST 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 so-called multi-conformer TST), which requires one to associate a unique reactant conformer with each corresponding transition state conformer, MS-VTST does not require this one-to-one 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 MS-VTST as the correct formulation of transition state theory with multiple conformers.

3.3 Multipath variational transition state theory

MS-VTST only considers one single lowest-energy 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. Multi-path variational transition state theory149,150 (MP-VTST) is a generalization based on MS-VTST.

We now derive the MP-VTST rate constant for a unimolecular reaction (for a bimolecular reaction the derivation follows a similar manner). For a unimolecular reaction, the MP-CVT 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:

 
image file: c7cs00602k-t104.tif(129)
where V1 is the classical barrier height of the lowest-energy reaction path; Q‡-SS-Tp is the single-structural rovibrational torsional partition function for the saddle point of path p, which represents the contribution from the p-th saddle point structure to the total conformational rovibrational torsional partition function of the transition state:
 
image file: c7cs00602k-t105.tif(130)
where Up 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 lowest-energy saddle point structure. Notice that the differences between the barrier heights of various paths have already been included in the Boltzmann factor of the single-structural rovibrational torsional partition function of the saddle point of path p in eqn (130), and therefore in eqn (129) only the V1 of the lowest-energy path was explicitly written out.

Eqn (129) can be further modified as follows

 
image file: c7cs00602k-t106.tif(131)
and it is equal to
 
image file: c7cs00602k-t107.tif(132)
where Q‡-SSHO1 is the single-structural rovibrational partition function of the saddle point of the lowest-energy path. Notice that if we included all the reaction paths, then the denominator of the first term in eqn (132) is equal to the MS-T conformational rovibrational partition function Q‡-MS-T. Finally, eqn (129) can be expressed in the following form
 
kMP-CVT/T(T) = 〈γFMS-TfwdkTST1(133)
where kTST1 is the conventional (single-structural) TST rate constant computed based on the lowest-energy path. The averaged generalized transmission coefficient 〈γ〉 is defined as
 
image file: c7cs00602k-t108.tif(134)
The multistructural and torsional anharmonicity factor FMS-Tfwd is
 
image file: c7cs00602k-t109.tif(135)
As we can see from eqn (133), if we only include the single lowest-energy path, the MP-VTST rate constant reduces to the MS-VTST 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 single-structural HO rovibrational partition functions QX-SSHO (X is saddle point, reactant, or variational transition state) in the above equations with the single-structural rovibrational partition function with torsional anharmonicity QX-SS-T;151 the computation of variational transmission coefficient is also based on QX-SS-T. For the generalized transition state, QGT-SS-T(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 ΓCVTp) will be different from the previous approach, because the generalized-transition-state Gibbs free energy of activation ΔGGT,0act(T,s) is now evaluated with the torsional anharmonicity; and hence the position of the generalized-transition-state dividing surface, where ΔGGT,0act(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 QX-SSHO. As a further complete treatment, in the full MP-VTST, we use QX-MS-T(s) to compute the variational transmission coefficient, which is computationally prohibitive.

The microcanonical variational transition state theory can also be generalized to include MS-T effects; this is multistructural microcanonical VTST (MS-μVT). And this is done by

 
image file: c7cs00602k-t110.tif(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 QMS-T.

3.4 Truncation of paths included in MP-VTST

Ideally, all the reaction paths generated by the multiple conformational structures of the transition state structure should be included in the computation of MP-VTST. However, for even a medium-sized 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 MP-VTST calculation, that is to say, in the practical MP-VTST calculation, P < N. If P = N, then the denominator of eqn (134) is exactly equal to QMS-T(TS); and if P < N, then we have introduced truncation errors in the computed MP-VTST rate constant.

A recent study152 compares the MP-VTST results including different numbers of paths to the MP-VTST with all the paths included; the comparison is for the hydrogen abstraction reaction of tert-butanol. 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 path-dependent, and higher-energy 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 VGa 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 higher-energy 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 MP-VTST is good enough (as compared to MP-VTST 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 one-dimensional 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 reaction-path VTST (RP-VTST) 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. Reaction-path VTST is not usually valid for a loose transition state because the accessible coordinate has wide-amplitude motion that cannot be well descried by torsions and nearly harmonic stretches and bends. It is still possible to use reaction-path 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 variable-reaction-coordinate VTST155–159 (VRC-VTST), as discussed next.

4.1 Variable-reaction-coordinate variational transition state theory

Next we review the VRC-VTST theory. The VRC-VTST methodology is available in the Polyrate code42 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 VRC-VTST, 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 high-pressure-limit rate constant given by VRC-VTST for association (R1 + R2 → P) is

 
image file: c7cs00602k-t111.tif(137)
where ge 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; Q1 and Q2 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 VRC-VTST is different from the one in reaction-path VTST. In VRC-VTST, the reaction coordinate is defined by pivot points. Pivot points are used to orientate the two fragments in the VRC-VTST 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 temperature-dependent. The reader is referred to the SI of our paper161 on 2CF2 ⇌ C2F4 for details of the way we use pivot points to define a reaction coordinate.

In the single faceted VRC-VTST, each reactant has one single pivot point, and we denote their coordinates as P1 and P2; the reaction coordinate s is defined as s = |P1P2|. 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 multi-faceted VRC-VTST, 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 rij, 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 rij set:

 
s = min{rij}(138)

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) = 〈Nq(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 2CF2 ⇌ C2F4.161 The number of states at a given orientation is computed as
 
image file: c7cs00602k-t112.tif(140)
 
image file: c7cs00602k-t113.tif(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 = v1 + v2 + 3, where v1 and v2 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 Ii,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 ka 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 VRC-VTST, the optimization of the transition-state 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 CVT-VRC-VTST, 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 μVT-VRC-VTST, 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 Nq(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 ensemble-averaged 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 long-range form of the interaction potential, although in general – using the language of association reactions – one must also consider the transformation of long-range 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 long-range 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 theory22,163 (CUS) has been developed as an extension to canonical ensembles of Miller's unified statistical theory164,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 re-dissociate to reactants. In the case of two transition states in series, the CUS rate constant is:

 
kCUS(T) = kCVT(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:
 
image file: c7cs00602k-t114.tif(143)
where qCVT(T) is the rovibrational partition function for the canonical variational transition state (the location of which is denoted as image file: c7cs00602k-t115.tif, which corresponds to the highest maximum point on the generalized-transition-state Gibbs free energy surface), qmax(T) is for the second highest maximum point on the generalized-transition-state Gibbs free energy surface (the location of which is denoted as image file: c7cs00602k-t116.tif), and qmin is for the lowest local minimum point that is between image file: c7cs00602k-t117.tif and image file: c7cs00602k-t118.tif on the generalized-transition-state Gibbs free energy surface. All these three partition functions are computed based on the zero-of-energy being at the minimum of reactants’ potential well (i.e., the value of VMEP at s = −∞); notice that, there is a relation between the generalized-transition-state vibrational partition function QGT(T,s) that is based on the zero of energy being at the value of VMEP(s) of the potential on MEP, and the vibrational partition function qGT(T,s) that is based on the former definition of zero of energy:
 
qGT(T,s) = exp[−VMEP(s)/kBT]QGT(T,s)(144)

Clearly, if image file: c7cs00602k-t119.tif, image file: c7cs00602k-t120.tif, and image file: c7cs00602k-t121.tif are the same, as in the case which there is only one local maximum on the generalized-transition-state 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 collision-induced 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 kCUS can also be expressed by using the rate constant (or the reactive flux coefficient) of passing the highest variational transition state kI, the second highest variational transition state kII, and the formation of the lowest local minimum complex kC, on the Gibbs free energy surface along the MEP, as follows:

 
image file: c7cs00602k-t122.tif(145)

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

 
image file: c7cs00602k-t123.tif(146)
where kIIi is the rate constant for the i-th 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 paths171–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 SN2 ion–molecule reactions,247 and for the competition between E2 and SN2 reactions.248 A non-statistical extension of CCUS model will be discussed in Section 6.

5. Pressure-dependent 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 high-pressure-limit rate constants. However, in practical laboratory measurements, the pressure is seldom high enough to measure the high-pressure limit, and measured rate constants can be significantly different from the ones computed for the high-pressure-limit. For thermal unimolecular reactions, the rate falls off as pressure decreases, and this is called the falloff effect. It is a consequence of a non-equilibrium state distribution in the reactant.

A master equation34,178–182 may be used to accurately treat non-equilibrium effects if the full matrix of state-selected reaction rates and state-to-state 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 state-specific 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 one-dimensional, and it may be written:183,184

 
image file: c7cs00602k-t124.tif(147)
 
image file: c7cs00602k-t125.tif(148)
where pi(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, ynA(t) and ynB(t) are concentrations of the nth arrangement-channel reactants with A and B distinguishing between the two reactants, kij(E), kin(E), and kni(E) are the microcanonical rate constants for isomerization, association and dissociation reactions, bn(E,T) is the Boltzmann distribution for bimolecular reactant channel n, and Nisom, Nreact, and Nprod are the numbers of isomers, bimolecular reactant channels, and bimolecular product channels, respectively.

Popular master equation solvers are available, such as Mesmer,185Multiwell,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 three-parameter Arrhenius expression by using Variflex, which includes variational effects, and then using the inverse Laplace transform method for the generalized 3-parameter 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.

Pressure-dependent 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 low-pressure rate constants of the formation of the stabilized adduct are lower than the high-pressure-limit, and the rate constants of the further isomerization/dissociation of the formed adduct are larger than the high-pressure equilibrium rate constant. Experimentally, Rabinovitch and co-workers 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 low-pressure rate constants are smaller than the high-pressure 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) theory189 or via the full microcanonical variational transition state theory (sometimes called variational RRKM). Barker, Westmoreland, Dean, and Bozzelli and their coworkers190–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 system-specific quantum RRK theory78,198,199 (SS-QRRK) 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 high-pressure-limit rate constants. In conventional RRKM theory for k(E), energy-resolved variational effects, multi-dimensional 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. SS-QRRK can be viewed as an effective way of recovering the above-mentioned effects in k(E) from the high-pressure-limit canonical k(T) that do include those effects in the canonical ensemble; and we use it to compute the pressure-dependence so that the computed k(T,p) can connect with the correct high-pressure-limit.

5.2 Thermal activation

The thermal activation mechanism for unimolecular reaction X → P is
 
image file: c7cs00602k-t126.tif(149)
 
image file: c7cs00602k-t127.tif(150)
where X is the reactant, M is bath gas (such as He, Ar, N2, 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 E0 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 kact, and we will call the reverse step de-activation and label it kdeact, 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 de-energizes X*. The rate constant for the step in which reactant is converted to product is called kcon,X (this could be either isomerization or dissociation). The phenomenological unimolecular reaction rate constant is defined as
 
image file: c7cs00602k-t128.tif(151)

In the classical Lindemann theory, all of the rate constants are considered to be energy-independent. If we use steady-state approximation to treat X*, i.e., if we assume

 
image file: c7cs00602k-t129.tif(152)
then we get
 
image file: c7cs00602k-t130.tif(153)

If we assume ideal-gas behavior, [M] can be replaced by p/RT, where p is pressure. Then, since kuni[X] = kcon,X[X*], we have

 
image file: c7cs00602k-t131.tif(154)

In the high-pressure limit, the Lindemann unimolecular rate constant becomes

 
kuni = kactkcon,X/kdeact(155)
and the unimolecular reaction is of pseudo-first order; whereas in the limit of very low pressure, the low-pressure-limit k0uni is much smaller than kuni:
 
k0uni = (p/RT)kact(156)
which depends linearly on the pressure, and the reaction is of second order. The falling off of the kuni 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 kcon,X should be considered. One could use microcanonical VTST to calculate the energy-dependent rate constants, just as RRKM theory uses microcanonical conventional TST energy-dependent rate constants. Alternatively one could calculate thermal rate constants and obtain the fixed-energy 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 SS-QRRK as an effective kinetics-specific inverse Laplace transform to convert canonical kcon,X rate constants to microcanonical ones; SS-QRRK 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 SS-QRRK, is that we neglect rotational energy and assume the vibrational energy E is randomly distributed among all the modes (nevertheless, in SS-QRRK the rotational contribution is implicitly included in k(E) because the canonical high-pressure-limit includes rotational modes, and thus they are implicitly in AK and E0, which are discussed below). For simplicity we treat X and X* as multi-mode 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 kinetics-specific 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 [small nu, Greek, macron] (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[v with combining macron](157)
where c is speed of light. Following a similar derivation to that in classical Lindemann theory, but employing n-dependent rate constants, the unimolecular rate constant is
 
image file: c7cs00602k-t132.tif(158)
in which m is the number of quanta at the threshold energy E0,
 
m = E0/hc[v with combining macron](159)
and the determination of threshold energy in SS-QRRK will be discussed later.

The ratio of the thermal activation and de-activation rate constants in eqn (149) is a mixed-ensemble 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 by200

 
image file: c7cs00602k-t133.tif(160)
where s is the number of vibrational degrees of freedom of X* (which is 3N − 6 for a nonlinear molecule). We assume that kdeact 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 kact, which yields
 
image file: c7cs00602k-t134.tif(161)

In QRRK, a unimolecular reaction happens when a specific vibrational mode associated with the reaction coordinate possesses threshold energy E0. Therefore the QRRK rate constant equals a frequency factor AK (the superscript stands for Kassel) times the fraction of molecules that have at least m quanta of excitation; this is given by181,201

 
image file: c7cs00602k-t135.tif(162)

Substituting eqn (159), (161) and (162) into (158) yields the QRRK rate constant, but in order to finish the calculation we must specify AK, E0, and kdeact. The computation of kdeact will be discussed in Section 5.4, and next we turn our attention to AK and E0.

Taking the high-pressure-limit of eqn (158), using eqn (159), and summing the series yields the following pure Arrhenius form:

 
kuni(T) = AK[thin space (1/6-em)]exp(−E0/kBT)(163)
and this motivates the approach we proposed to obtain the parameters needed in SS-QRRK computations. The first step is to calculate the high-pressure-limit canonical-ensemble VTST rate constant (for example, an MP-CVT/SCT calculation that includes variational effects, multidimensional tunneling, various anharmonicity). Next, we fit the canonical-ensemble high-pressure-limit rate constant by the four-parameter formula202
 
image file: c7cs00602k-t136.tif(164)
where B, T0, E′, and n′ are fitting constants. In combustion community, the widely used fitting formulas are either conventional Arrhenius formula or the three-parameter fitting formula in the form BTn[thin space (1/6-em)]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 non-linearly reduce the activation energy significantly,150 and also at high temperature. Our four-parameter formulas are more physically sound, and they permits a nonzero rate constant for an exothermic reaction in the low-temperature limit; furthermore they are more accurate in terms of the fitting error.202

Based on our four-parameter 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,

 
image file: c7cs00602k-t137.tif(165)
which yields
 
image file: c7cs00602k-t138.tif(166)

With this energy of activation, we have a local Arrhenius fit for temperature T given by

k(T) = A(T)exp[−Ea(T)/RT]
with A(T) obtained by
 
A(T) = k(T)exp[Ea(T)/RT](167)

Then Ea(T) and A(T) are used for E0 and AK in the QRRK expression given by eqn (158), (159), (161) and (162). This yields SS-QRRK. Notice that the parameters E0 and AK are constants in conventional QRRK, but they depend on temperature in SS-QRRK, 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 high-pressure-limit rate constants and approximately transfer it to the microcanonical energy-dependent rate constants, without actually doing the microcanonical computations. We view the subsequently obtained SS-QRRK 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 MS-T anharmonicity and microcanonical multidimensional tunneling), the sum in eqn (158) can be replaced by an integration, which yields

 
image file: c7cs00602k-t139.tif(168)
in which
 
image file: c7cs00602k-t140.tif(169)
where ρ(E) is MS-T rovibrational density of states of X, QX(T) is the MS-T rovibrational partition function of reactant X, Ethr 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 ZPE-included electronic potential-energy barrier height; this is different from the way of determining the threshold energy in SS-QRRK theory, which is a temperature-dependent 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 small-curvature tunneling (SCT), the MS-μVT/SCT rate constant is
 
image file: c7cs00602k-t141.tif(170)
where the cumulative reaction probability N‡,MS-T/SCT is the convolution of the SCT tunneling probability PSCT(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 ground-state 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 medium-sized systems.

Finally, let us consider use SS-QRRK for treating unimolecular dissociation with multiple parallel dissociation reactions.

 
image file: c7cs00602k-t142.tif(171)
 
image file: c7cs00602k-t143.tif(172)
where the reactant X has multiple parallel dissociation reactions, of which the products are respectively P1, P2,…, PN (notice that here the notation Pi for the i-th 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:
 
image file: c7cs00602k-t144.tif(173)
which, by using the condition [X] = [P1] + [P2] + ⋯ + [PN], can be re-written as the summation of the rate constants of all the parallel dissociation reactions
 
image file: c7cs00602k-t145.tif(174)
where
 
image file: c7cs00602k-t146.tif(175)

And the SS-QRRK pressure-dependent rate constant for the i-th dissociation reaction, which is a straightforward generalization of eqn (158), is:

 
image file: c7cs00602k-t147.tif(176)
where mi 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 fall-off in the cases such as SiH4 → SiH3 + H or SiH2 + H2, although theirs work was based on conventional TST without including the dynamic effects we consider in SS-QRRK.

The addition of rate constants for parallel reactions, as in eqn (174), is clearly valid for microcanonical rate constants and for high-pressure-limit 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 fall-off region if high-energy states are depleted more rapidly by reactions with lower-energy thresholds than they are by reactions with high-energy thresholds. This phenomenon has been discussed in conjunction with the treatment of the unimolecular reactions of 2-methylhexyl 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 two-step reaction with a unimolecular intermediate, the chemical activation mechanisms are:

(a) association reaction (Y + Z → YZ)

 
image file: c7cs00602k-t148.tif(177)

(b) two-step reaction (Y + Z → YZ → P)

 
image file: c7cs00602k-t149.tif(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 two-step reaction), P is the product of further reaction of the intermediate YZ*. In chemical activation, the rate constant kadd of the formation of YZ* is not pressure-dependent because it is formed via a bimolecular chemical reaction, of which the value is just the high-pressure-limit temperature-dependent bimolecular reaction rate kadd(T). The rate constant for the dissociation of YZ* back to the reactants Y and Z is krev, which is treated as energy-dependent; the rate constant for the formation of the product(s) from YZ* is kcon,YZ(E).

Following a similar treatment as for thermally activated unimolecular reactions, and defining kstab as the phenomenological bimolecular stabilization rate constant,

 
image file: c7cs00602k-t150.tif(179)
and krxn as the phenomenological bimolecular rate constant for reaction
 
image file: c7cs00602k-t151.tif(180)

We applied steady-state approximation to the activated adduct YZ*. Then the pressure-dependent rate constants given by SS-QRRK theory are:

(a) for association reaction

 
image file: c7cs00602k-t152.tif(181)

(b) for two-step reaction

 
image file: c7cs00602k-t153.tif(182)
and
 
image file: c7cs00602k-t154.tif(183)
where f(E) is the fraction of energized adduct (YZ*) at energy E, which is given by
 
image file: c7cs00602k-t155.tif(184)
and Kact,YZ is given by eqn (157) and (160), and E0 is the SS-QRRK threshold energy for YZ dissociating back to reactants Y and Z, and krev and kcon,YZ are computed by eqn (162) with the appropriate m (number of quanta at threshold energy) and AK values (for krev, these parameters are computed from the reaction of YZ dissociating to Y and Z; and for kcon,YZ, from the reaction of YZ dissociation to P).

Instead of using SS-QRRK 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

 
image file: c7cs00602k-t156.tif(185)

(b) two-step reaction

 
image file: c7cs00602k-t157.tif(186)
 
image file: c7cs00602k-t158.tif(187)
where the energy distribution function
 
image file: c7cs00602k-t159.tif(188)
and in μVT, the threshold energy is denoted as Ethr, and it is ZPE-included electronic potential-energy barrier height, which is different from the threshold energy used in SS-QRRK theory (which is a temperature-dependent 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 SS-QRRK for chemical activation processes with the intermediate YZ undergoing multiple further parallel dissociation reactions. The mechanism is as follows:

 
image file: c7cs00602k-t160.tif(189)
 
image file: c7cs00602k-t161.tif(190)

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

 
image file: c7cs00602k-t162.tif(191)

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

 
image file: c7cs00602k-t163.tif(192)
where E0 is the SS-QRRK 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 high-energy states by the low-threshold and high-threshold reactions are comparable, then the current treatment for parallel channels provides a good approximation. And the pressure-dependent rate constant k(i)rxn for the formation of the product(s) Pi is:
 
image file: c7cs00602k-t164.tif(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., de-activation has occurred). Note that de-activation refers to reducing the total energy below the threshold for reaction. However, in reality it may takes multiple energy transfer collisions to fully de-activate 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 single-well, single-reaction-channel systems by a master-equation-based weak-collision 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 kdeact (in units of cm3 molecule−1 s−1) is computed as:
 
kdeact = βcΩ2,2*kHS(194)
where kHS is the hard-sphere 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
 
image file: c7cs00602k-t165.tif(195)
where εA–M is the Lennard-Jones 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:
 
image file: c7cs00602k-t166.tif(196)
where |〈ΔE〉| is the value of the average energy transferred during both energization and de-energization processes (this quantity is also denoted as |〈ΔEall|, in order to be distinguishable from the energy transferred only for de-energization 〈ΔEdown,179,210,211 which is a positive number) and
 
image file: c7cs00602k-t167.tif(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 〈ΔEdown parameter is available, then the collision efficiency could be computed as:179
 
image file: c7cs00602k-t168.tif(198)

The FE0 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 FE0 by using the MS-T partition function, which is reviewed in Section 3.1, in the inverse Laplace transform; if the single structure quasi-harmonic oscillator model is a good approximation, FE0 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 |〈ΔEall| 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 so-obtained pressure-dependence model becomes predictive for other temperatures and pressures.

For the numerical value of FE0, 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, FE0 can deviate from 1 or 1.15 by quite a large amount; in fact, for a large-sized molecule with lots of vibrational degrees of freedom at high temperatures, FE0 grows exponentially. For instance, for C2F4 dissociating to two CF2,161FE0 is computed to be 1.01 at 300 K, 1.12 at 400 K, and 1.39 at 1000 K; for (SiH3)2SiHSiH isomerization to (SiH3)2SiSiH3,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 done189,190 in the literature. When the density of states is obtained via our MS-T partition function (see Section 3.1) and FE0 is computed by eqn (197), the multiple structural and torsional anharmonicity (and also, the vibrational anharmonicity for high-frequency modes when a vibrational scale factor is applied) is built into FE0, and therefore the previously proposed empirical internal rotor correction and anharmonicity correction212 to FE0 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 |〈ΔEall|, 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, |〈ΔEall| 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 literature222–228 for examples of experimental values of |〈ΔEall|.

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 SS-QRRK calculations for the collision efficiency βc. For H addition on toluene (Ar as bath gas, with the energy transfer parameter |〈ΔEall| 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 SO2 and OH association reaction in the atmosphere (assuming N2 as bath gas, with the energy transfer parameter |〈ΔEall| being 74 cm−1),232 it is 0.21 at 250 K, 0.18 at 298 K, and 0.14 at 400 K. For C2F4 dissociating to two CF2 (with Ar as bath gas, with the energy transfer parameter |〈ΔEall| 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 pressure-dependent rate constant k(T,p)

Functional forms have previously been proposed for fitting the pressure-dependent 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 one-dimensional fits as a function of temperature.

Step 1: fit the computed high-pressure-limit rate constant k(T) (for a bimolecular association reaction, the units are cm3 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 low-pressure-limit rate constant k0(T), which is defined as the value of k(T,p)/[M] in the limit of p going to zero. The low-pressure-limit rate constant is a pseudo-third-order rate constant with unit being cm6 molecule−2 s−1 for a bimolecular association, and it is a pseudo-second-order rate constant with unit being cm3 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.,

 
image file: c7cs00602k-t169.tif(199)
where A, n, E and T0 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 p1/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 high-pressure-limit. The p1/2(T) could be fit with many possible formulas, and we give one of the fitting formulas here:

 
image file: c7cs00602k-t170.tif(200)
where a1, a2, a3, l1, l2, T1, and T2 are the fitting parameters, in which a1, a2, and a3 are unitless, and l1, l2, T1, and T2 are in K.

The final fitted pressure-dependent rate constant k(T,p) is expressed in terms of the fitted k(T), k0(T), and p1/2(T) (and thus it does not require any additional fitting) by the following interpolatory equation:

 
image file: c7cs00602k-t171.tif(201)
where
 
image file: c7cs00602k-t172.tif(202)
and
 
image file: c7cs00602k-t173.tif(203)
in which the Boltzmann constant kB in eqn (202) is 1.38 × 10−22 bar cm3 molecule−1 K−1, p is pressure in bar, B is a derived parameter in units of bar−1, and d is in bar2. Eqn (201) has the following three properties, which ensure the fitted rate constants have the right high-pressure-limit, low-pressure-limit, and the transition pressures:
 
k(T,p = ∞) = k(T)(204)
 
k(T,p = 0) = k0(T)[M](205)
and
 
k[T,p = p1/2(T)] = k(T)/2(206)

The pressure-dependent activation energy Ea(T,p) at a given pressure p = p0 can then be determined by numerical differentiation by employing the following formula:

 
image file: c7cs00602k-t174.tif(207)

In the high-pressure 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 temperature-dependent rate constants as in a canonical ensemble; microcanonical equilibrium when we discuss energy-dependent 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 local-equilibrium or quasi-equilibrium assumption and the no-recrossing assumption. We have stated the assumptions in classical mechanical language, but they can be rephrased quantum mechanically by using Heisenberg-transformed 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 mechanics234 (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 kinetic-isotope-effects literature,235–237 but we prefer to save the word semiclassical for WKB-like treatments of tunneling. Thus when we say “semiclassical,” we refer to including tunneling, but when they say “semiclassical” in the kinetic-isotope-effects 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 (quasi-equilibrium, 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 + H2 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 three-dimensional 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 + CH4 → H2 + CH3 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 high-pressure-limit rate constant of an ensemble, where the high-pressure limit is invoked to main the local equilibrium of reactant states. Strictly speaking, non-equilibrium effects are not included, although sometimes one extends the theory by including a non-equilibrium 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 pressure-dependent 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 local-equilibrium 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 low-temperature 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 quasi-equilibrium 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 vibrational-energy relaxation (“fast IVR”), and fast IVR may be considered part of the transition state theory quasi-equilibrium assumption in this context. Thus transition state theory, and hence also RRKM theory, is not applicable to very short-lived complexes that do not have time for full intramolecular energy redistribution241 or during the induction period of a long-lived complex, i.e., during that period before a steady-state phenomenological rate constant is established.

Next we turn to the no-recrossing assumption.

Hase and co-workers 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 transition-state theory may be quite different in classical and quantum mechanics.26 Interestingly, a canonical unified statistical theory163 (which is an extension of the earlier unified statistical theory245,246) has been applied to study similar problems,247,248 and the computed kinetic isotopic effects (KIEs) agree reasonably well with experimental observations for gas-phase SN2 reactions.

Consider, as an example, association reactions involving small molecules, where the system crosses a loose bottleneck. If there are not enough low-frequency 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 small-molecule 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 high-free-energy 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 beyond-transition-state-theory assumptions. For example, we combined VTST with non-statistical assumptions to study bifurcation in the BH3 + propene system, and we found that our computed selectivity agrees well with experimental observations.257 This method is called the canonical competitive non-statistical (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 non-statistical 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 ktot, which is the sum of the direct (kdir) and indirect (kind) components, is computed as:

 
image file: c7cs00602k-t175.tif(208)
where k*2 = min(kC,k*2A + k*2B), and ki 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 (kdir) is computed as:
 
kdir = k*1k*2/kC(209)
and the total reaction rate constants for the two branches A and B are:
 
image file: c7cs00602k-t176.tif(210)
 
image file: c7cs00602k-t177.tif(211)
where RA/B is the branching ratio of the direct reactive components, which is approximated as the equilibrium constant KA/B given by statistical phase space theory260,261 with additional non-statistical corrections:
 
image file: c7cs00602k-t178.tif(212)
where the rovibrational partition function QX (X = A or B, which denote the final product of reaction A or reaction B) is
 
image file: c7cs00602k-t179.tif(213)
in which QX,⊥ is the rovibrational partition function that includes all the bound-mode vibrations that are orthogonal to the reaction coordinate s, but excluding the reaction coordinate s itself, which becomes a bound-mode vibration in the product (and it is denoted as mode F, and its corresponding vibrational quantum number is nF 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 available258 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 (sproductsTS1) in coordinates scaled to a reduced mass to a characteristic velocity in the reaction coordinate for reaction coordinate quantum number nF.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 high-level 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 co-workers262 presented an interesting study for the same BH3 + 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 temperature-dependence of gas-phase 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 〈ΔEdown are 900–1100 cm−1. We may compare this to the gas phase where the usually used value of 〈ΔEdown 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 pressure-dependent VTST calculations with SS-QRRK theory. Rigorously speaking, the pressure-dependence calculation is beyond transition state theory; transition state theory just provides the key ingredients (high-pressure-limit rate constant, and the effectively transformed energy-resolved microcanonical SS-QRRK rate constants) that are needed in the computation of pressure-dependent rate constant.

One important assumption about the non-statistical 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 exponential-down 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 CHF3 dissociation reaction198 (the dissociating products are HF and CF2) and the C2F4 dissociation reaction (which is a uphill from reactant to two CF2 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 SS-QRRK 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 J-dependence in such a SS-QRRK 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) SS-QRRK is an efficient model to approximately capture all the physical effects contained in high-pressure-limit 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 ensemble-averaged 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 so-called mixed discrete-continuum 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

 
image file: c7cs00602k-t180.tif(214)
where the volume integral is carried out over the range of R corresponding to species S; and C is a geometry-independent constant. We can carry out this integral in two steps: first
 
image file: c7cs00602k-t181.tif(215)
and then
 
image file: c7cs00602k-t182.tif(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
 
image file: c7cs00602k-t183.tif(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) + ΔGS*(R)(218)
where R denotes the 3N − 6 internal solute coordinate, V(R) is the gas-phase Born–Oppenheimer electronic potential energy surface, and ΔGS*(R) is the fixed-concentration 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 liquid-phase solution; that is an ideal gas at a gas-phase 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 SMD278), ΔGS* has already been included in the SCF single-point energy, i.e., what we have calculated directly from a self-consistent reaction field (SCRF) calculation is the free energy surface.

The molar Gibbs free energy of a species i in liquid-phase solution is:277

 
Gi = Gi* − TSlib,i(219)
where Gi* 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 pseudo-chemical potential, and it is associated with the internal free energy of the species in solution including its coupling to the solvent); one may approximate Gi* in as:
 
Gi* = U0 + Gint(220)
where Gint is the internal thermal Gibbs free energy summed over all the configurations; U0 is the zero-point inclusive free energy surface, given in the harmonic approximation by:
 
image file: c7cs00602k-t184.tif(221)
where Fvib is the number of vibrational degrees of freedom, and ωm is the vibrational frequency of mode m. The entropy term Slib,i in eqn (219) is the entropy of liberation of species i, which is given by:
 
image file: c7cs00602k-t185.tif(222)
where R is the gas constant, ci is concentration of i in moles per unit volume, NA 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 standard-state Gibbs free energy (which is directly related to partition functions and standard equilibrium constants), the standard-state Gibbs free energy of species B in solution image file: c7cs00602k-t186.tif (i.e., the chemical potential of B) should be computed as:

 
image file: c7cs00602k-t187.tif(223)
where image file: c7cs00602k-t188.tif is the standard-state gas-phase Gibbs free energy of species B, which is usually computed based on gas-phase computed frequencies of the gas-phase optimized geometries; but one can also use solution-phase computed frequencies of solution-phase optimized geometries.279 The superscript “°” in eqn (223) denotes the usual thermodynamic standard state, which is defined as follows: (a) gas-phase 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 liquid-phase 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 image file: c7cs00602k-t189.tif in eqn (223) is the amount that one must add to the fixed-concentration 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 right-hand side of the solvation process (the liquid-phase solution) at 1 mol L−1; the value of image file: c7cs00602k-t190.tif 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 condensed-phase 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 non-negligible. 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 liquid-phase reaction:281,289 (1) separable equilibrium solvation (SES); (2) equilibrium solvation path (ESP); (3) non-equilibrium 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 gas-phase, and then adds the free energy of solvation (at a fixed solute geometry) to the gas-phase free energy for each point on the gas-phase reaction path, using the gas-phase optimized geometries and gas-phase 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 gas-phase 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 ESP281 (within a factor of 2, when the same solvation model269,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.

Non-equilibrium solvation can be considered via a linear-response theory for solvent friction effects, which can be incorporated with VTST and tunneling;290,291 the key concept in such linear-response 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 Doll292 compared classical trajectory simulations results with classical TST results for diffusion constants. In trajectory simulations, the diffusion constants D are obtained by:
 
image file: c7cs00602k-t191.tif(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 〈Δr2(t)〉 is the mean-squared displacement of the migrating atom, which can be obtained by averaging over trajectories. In practice, D is obtained by plotting 〈Δr2(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 Arrhenius-type equation:
 
D = D0eEa/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

 
image file: c7cs00602k-t192.tif(226)
where l is the distance (hop length) between adjacent minimum-energy binding sites, and khop is the hopping rate. One-dimensional harmonic conventional classical transition state theory yields
 
khop = npv0eV/kBT(227)
where V is the energy difference between saddle point (between two binding sites) and local binding site minimum; v0 is the harmonic frequency of the migrating atom moving along a pseudo-one-dimensional potential, in which all the other atoms are adiabatically relaxed; and np is the number of binding sites accessible by a single hop. If we equate V to Ea and equate eqn (225) and (226), the parameter D0 in the Arrhenius-type equation for the diffusion constant is:
 
image file: c7cs00602k-t193.tif(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:

 
image file: c7cs00602k-t194.tif(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 as295
 
δ[f(R)]|∇f|〉A = fB/ω(230)
where ω is the width of the simulation box, and fB 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 khop, and therefore it can give a better prediction for diffusion constants. The practical computations can be carried out with an embedded-cluster model.296–300 In such model systems, there is only one non-metal 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 non-fixed 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 next-door 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 observed305–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 temperature-independent due to tunneling. The quantitative definition of the transition temperature is the temperature at which the Arrhenius plot has maximum curvature. This zero-temperature-limit, 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. Reanalysis300 of the low-temperature situation shows that the experimental “transition temperature” actually signifies the transition from ground-state-dominated tunneling to thermally activated tunneling. Thermally activated tunneling continues to dominate the rate constant at temperatures far above the transition temperature. In low-temperature-limit, 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 harmonic-oscillator approximation, the low-temperature-limit diffusion coefficients are:

 
image file: c7cs00602k-t195.tif(231)
where l is the lateral distance between the two minimum-energy 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, [small upsilon, Greek, macron]R is the vibrational frequency in wavenumber for the lowest-frequency hydrogenic vibration at the reactant well (which corresponds to the reaction coordinate over most of the reaction path away from the reactant well), and PG(ER0) is the quantum transmission probability at the harmonic vibrational ground-state 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. Ensemble-averaged VTST309 (EA-VTST) with multidimensional tunneling was developed for modeling kinetics in enzymes, and we will briefly review EA-VTST 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

EA-VTST 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 pre-defined 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 WC(z) (free energy profile, or – more technically – the generalized free energy of activation profile). Once the classical potential of mean force is available, the generalized-transition-state (GT) classical potential of mean force activation energy ΔWGT,C(T,z = z*) at temperature T is computed by:
 
image file: c7cs00602k-t196.tif(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, ΔWGT,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 reaction-coordinate-mode free energies of the reactant to the computed classical ΔWGT,C(T,z = z*); and the so-obtained quantity is the generalized transition state classical free energy of activation given by

 
ΔGGT,C(T,z = z*) = ΔWGT,C(T,z = z*) − GC,F(T,z = zR)(233)
where GC,F(T,z = zR) is the classical free energy contribution of the reaction coordinate at its value at the reactant state zR; 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 N1-atom subsystem; in particular, only 3N1 − 7 modes are quantized at the transition state, and 3N1 − 6 modes are quantized at the reactant. The Gibbs free energy of activation obtained this way is called the single-reaction-coordinate (SRC) quasiclassical (QC) generalized-transition-state (GT) free energy of activation, which is denoted as ΔGSRCGT,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:

 
image file: c7cs00602k-t197.tif(234)

Stage 2. In this stage, we calculate ensemble-averaged recrossing and tunneling transmission coefficients in order to obtain the ensemble-averaged quasiclassical variational transition state theory rate constant. For EA-VTST with tunneling (abbreviated as “T”), the rate constant is:
 
kEA-VTST/T = γk(1)(235)
where k(1) is the VTST rate constant obtained in stage 1, and γ is ensemble-averaged transmission coefficient
 
image file: c7cs00602k-t198.tif(236)
where M is the total number of ensemble members included in the average, Γi is the recrossing transmission coefficient of member i, and κTi 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 secondary-zone approximation. This stage involves semiclassical reaction-path calculations for a subsystem (which is called the primary zone) within the secondary-zone 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 secondary-zone (also called “the frozen bath”) assumption works very well, and the non-equilibrium solvation effects are negligible; it is sufficient to stop at stage 2. As a further improvement, the entropic contributions from the secondary-zone can be included in the transmission coefficient. Then in the equilibrium secondary-zone approximation (which can be viewed as an analogous to a liquid-phase 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 secondary-zone 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 next-generation combustion engines. It has been demonstrated that anharmonicity150,327–330 plays an important role in controlling the VTST rate constants of hydrogen abstraction reactions. One interesting study that uses MS and MP-VTST in analyzing the effect of hydrogen-bonded transition states on rate constants (for gas-phase hydrogen abstraction reactions) shows that,150 the conventional thinking of H-bonded 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 HO2,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 radicals146,371–377 that are abundant in combustion flames. Guan and co-workers have also applied MS-VTST 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 engine-relevant conditions, have also been studied using VTST.380 Klippenstein and co-workers 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 pressure-dependence 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 climate-controlling atmospheric aerosol formation.

We have modeled the dissociation of C2F4 with VRC-VTST and SS-QRRK 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 SO2 with OH in the atmosphere;232 this reaction plays an essential role in acid rain and it is still attractive to experimentalists up-to-date,386 and our temperature- and pressure-dependent rate constant computations help to resolve the discrepancies between various experimental studies. Cartoni et al. have studied the formation of HSO2+ ion in the gas phase, which is of special interest because of an interesting reverse temperature-dependent kinetics.387 In the study of the important roles played by water and ammonia in promoting atmospheric reaction in atmosphere: the rate constant of water-catalyzed H2SO4 + 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 + H2S reaction has been understood, and predictions have been made for this reaction at higher temperatures.390 Loerting and co-workers have studied the proton transfer in small cyclic water clusters, which are important in understanding water-cluster 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 co-workers400 have studied the heterogeneous reactions for chlorine nitrate hydrolysis on water-ice 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

Silicon-based 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 MS-T 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 group-additivity 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 silicon-based 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 silane-based plasmas, and good agreements with experimental values was achieved.

8.4 Organic chemistry in solution

VTST with state-of-the-art solvation models can be used to calculate reaction rates and elucidate reaction mechanisms in liquid-phase 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 second-order mechanism for β eliminations and non-concerted mechanisms with discrete carbanion intermediates; distinguishing these experimentally is very difficult.

Liquid-phase 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 functional415 yields more accurate rate constants than does M08-HX, which is a density functional that is generally more accurate for gas-phase reactions; this is possibly because the error introduced in the solvation model partially cancels the error in the less-accurate 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

Ensemble-averaged 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 short-chain acyl-CoA dehydrogenase,427 catechol o-methyltransferase,428 glyoxalase I,429 coenzyme B12,430,431 HIV-1 protease,432 4-oxalocrotonate tautomerase enzyme,433 hydride transfer between Anabaena Tyr303Ser FNRrd/FNRox 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 tunneling239,448 (as opposed to one-dimensional 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 multiple-structure and torsional-potential anharmonicity in VTST. We reviewed variable-reaction-coordinate VTST as well as reaction-path VTST. We have discussed the recent extension of VTST for the convenient calculation of pressure-dependent 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 reaction-path variational transition state theory (RP-VTST) with multidimensional tunneling (MT); for barrierless reactions, use variable reaction coordinate variational transition state theory (VRC-VTST). Carry out direct dynamics rate constant calculations using one of these methods.

Gas-phase reactions Liquid-phase reactions
3-1 For RP-VTST, if needed, include the effects of multistructural torsional anharmonicity (MS-T) in the partition functions; this may involve extensive conformational searches. Compute MS-VTST/MT rate constants. If affordable, compute MP-VTST/MT rate constants. For reactions in the condensed phase, employ the separable equilibrium solvation (SES) VTST, the equilibrium solvation path (ESP) VTST, or the ensemble-average (EA) VTST.

For enzyme reactions, employ ensemble-averaged VTST.

3-2 For reactions whose rate constants are pressure-dependent, do pressure-dependence calculations either by the SS-QRRK method based on canonical VTST rate constants or by using a master equation solver with microcanonical VTST rate constants.
4 Compare the computed rate constants with experiments and explain any differences.
5 For the purpose of studying a complex reaction network, add the computed elementary reaction rate constants to the reaction model or update the older reaction kinetics data in the model.

Variational transition state theory, combined with the power of modern electronic structure methods, has been applied successfully to study the kinetics of a variety of chemical systems. We believe that the recent extensions of VTST can be useful tools for future research in combustion, atmospheric chemistry, environmental science, and other fields in which accurate kinetics data (which could be experimentally very hard to measure) play an indispensible role in understanding the detailed mechanisms of the chemical processes. VTST has also been applied successfully to enzyme kinetics and other condensed-phase reactions. With the recent progress of the development of more accurate and affordable electronic structure methods for treating inherently multiconfigurational systems, VTST is readily applicable to study the detailed kinetics in modern organometallic and inorganic catalysis.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported in part by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0015997. J. L. Bao also acknowledges a Doctoral Dissertation Fellowship (DDF) provided by University of Minnesota.

References

  1. P. Waage and C. M. Guldberg, Forhandlinger: Videnskabs-Selskabet i Christiania, 1864, vol. 35 (see translation in J. Chem. Educ., 1986, 63, 1044) Search PubMed.
  2. S. A. Arrhenius, Z. Phys. Chem., 1889, 4, 96 Search PubMed.
  3. S. Glasstone, K. Laidler and H. Eyring, The Theory of Rate Processes, McGraw-Hill, New York, 1941 Search PubMed.
  4. J. C. Polanyi and A. H. Zewail, Acc. Chem. Res., 1995, 28, 119 CrossRef CAS.
  5. R. Daudel, G. Leroy, D. Peeters and M. Sana, Quantum Chemistry, John Wiley & Sons, Chichester, 1983 Search PubMed.
  6. T. C. Allison and D. G. Truhlar, Testing the Accuracy of Practical Semiclassical Methods: Variational Transition State Theory with Optimized Multidimensional Tunneling, in Modern Methods for Multidimensional Dynamics Computations in Chemistry, ed. D. L. Thompson, World Scientific, Singapore, 1998, p. 618 Search PubMed.
  7. D. C. Chatfield, R. S. Friedman, D. G. Truhlar, B. C. Garrett and D. W. Schwenke, J. Am. Chem. Soc., 1991, 113, 486 CrossRef CAS.
  8. D. C. Chatfield, R. S. Friedman, D. W. Schwenke and D. G. Truhlar, J. Phys. Chem., 1992, 96, 2414 CrossRef CAS.
  9. D. C. Chatfield, R. S. Friedman, S. L. Mielke, G. C. Lynch, T. C. Allison, D. G. Truhlar and D. W. Schwenke, Computational Spectroscopy of the Transition State, in Dynamics of Molecules and Chemical Reactions, ed. R. E. Wyatt and J. Z. H. Zhang, Marcel Dekker, New York, 1996, p. 323 Search PubMed.
  10. D. C. Chatfield, R. S. Friedman, D. G. Truhlar and D. W. Schwenke, Faraday Discuss. Chem. Soc., 1991, 91, 289 RSC.
  11. R. S. Friedman and D. G. Truhlar, Chem. Phys. Lett., 1991, 183, 539 CrossRef CAS.
  12. R. S. Friedman, V. D. Hullinger and D. G. Truhlar, J. Phys. Chem., 1995, 99, 3184 CrossRef CAS.
  13. R. S. Friedman and D. G. Truhlar, Barrier Resonances and Chemical Reactivity, in Multiparticle Quantum Scattering with Applications to Nuclear, Atomic, and Molecular Physics, ed. D. G. Truhlar and B. Simon, Springer, New York, 1997, p. 243 Search PubMed.
  14. E. Wigner, J. Chem. Phys., 1937, 5, 720 CrossRef CAS.
  15. J. Horiuti, Bull. Chem. Soc. Jpn., 1938, 13, 210 CrossRef.
  16. J. C. Keck, J. Chem. Phys., 1960, 32, 1035 CrossRef CAS.
  17. J. C. Keck, Adv. Chem. Phys., 1967, 13, 85 CrossRef.
  18. B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 1979, 70, 1593 CrossRef CAS.
  19. B. C. Garrett and D. G. Truhlar, J. Phys. Chem., 1979, 83, 1052 CrossRef CAS.
  20. D. G. Truhlar and B. C. Garrett, Acc. Chem. Res., 1980, 13, 440 CrossRef CAS.
  21. D. G. Truhlar and B. C. Garrett, Annu. Rev. Phys. Chem., 1984, 35, 159 CrossRef CAS.
  22. D. G. Truhlar, A. D. Isaacson and B. C. Garrett, Generalized Transition State Theory, in Theory of Chemical Reaction Dynamics, ed. M. Baer, CRC Press, Boca Raton, 1985, p. 65 Search PubMed.
  23. B. C. Garrett and D. G. Truhlar, Transition State Theory, in Encyclopedia of Computational Chemistry, ed. P. V. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman and H. F. Schaefer III, John Wiley & Sons, Chichester, UK, 1998, vol. 5, p. 3094 Search PubMed.
  24. B. C. Garrett and D. G. Truhlar, Variational Transition State Theory, in Theory and Applications of Computational Chemistry: The First Forty Years, ed. C. E. Dykstra, G. Frenking, K. S. Kim and G. E. Scuseria, Elsevier, Amsterdam, 2005, p. 67 Search PubMed.
  25. A. Fernandez-Ramos, B. A. Ellingson, B. C. Garrett and D. G. Truhlar, Rev. Comput. Chem., 2007, 23, 125 CAS.
  26. D. G. Truhlar, W. L. Hase and J. T. Hynes, J. Phys. Chem., 1983, 87, 2664 CrossRef CAS.
  27. D. G. Truhlar, D.-h. Lu, S. C. Tucker, X. G. Zhao, A. Gonzàlez-Lafont, T. N. Truong, D. Maurice, Y.-P. Liu and G. C. Lynch, ACS Symp. Ser., 1992, 502, 16 CrossRef CAS.
  28. D. G. Truhlar, B. C. Garrett and S. J. Klippenstein, J. Phys. Chem., 1996, 100, 12771 CrossRef CAS.
  29. M. B. Hall, P. Margl, G. Náray-Szabó, V. L. Schramm, D. G. Truhlar, R. A. van Santen, A. Warshel and J. L. Whitten, ACS Symp. Ser., 1992, 721, 2 CrossRef.
  30. J. Gao and D. G. Truhlar, Annu. Rev. Phys. Chem., 2002, 53, 467 CrossRef CAS PubMed.
  31. D. G. Truhlar, J. Gao, M. Garcia-Viloca, C. Alhambra, J. Corchado, M. L. Sanchez and T. D. Poulsen, Int. J. Quantum Chem., 2004, 100, 1136 CrossRef CAS.
  32. M. Garcia-Viloca, J. Gao, M. Karplus and D. G. Truhlar, Science, 2004, 303, 186 CrossRef CAS PubMed.
  33. D. G. Truhlar, Variational Transition State Theory and Multidimensional Tunneling for Simple and Complex Reactions in the Gas Phase, Solids, Liquids, and Enzymes, in Isotope Effects in Chemistry and Biology, ed. A. Kohen and H.-H. Limbach, Marcel Dekker, Inc., New York, 2006, p. 579 Search PubMed.
  34. A. Fernández-Ramos, J. A. Miller, S. J. Klippenstein and D. G. Truhlar, Chem. Rev., 2006, 106, 4518 CrossRef PubMed.
  35. D. G. Truhlar and B. C. Garrett, Variational Transition State Theory, in The Treatment of Hydrogen Transfer Reactions, in HydrogenTransfer Reactions, ed. J. T. Hynes, J. P. Klinman, H. H. Limbach and R. L. Schowen, Wiley-VCH, Weinheim, Germany, 2007, vol. 2, p. 833 Search PubMed.
  36. A. Dybala-Defratyka, P. Paneth and D. G. Truhlar, Quantum Catalysis in Enzymes, in Quantum Tunneling in Enzyme-Catalyzed Reactions, ed. R. K. Allemann and N. S. Scrutton, RSC Publishing, Cambridge, UK, 2009, p. 36 Search PubMed.
  37. B. C. Garrett, D. G. Truhlar and R. S. Grev, Determination of the Bottleneck Regions of Potential Energy Surfaces for Atom Transfer Reactions by Variational Transition State Theory, in Potential Energy Surfaces and Dynamics Calculations, ed. D. G. Truhlar, Plenum Press, New York, 1981, p. 587 Search PubMed.
  38. D. G. Truhlar, A. D. Isaacson, R. T. Skodje and B. C. Garrett, J. Phys. Chem., 1982, 86, 2252 CrossRef CAS ; erratum: 1983, 87, 4554.
  39. M. M. Kreevoy and D. G. Truhlar, Transition State Theory, in Investigation of Rates and Mechanisms of Reactions, ed. C. F. Bernasconi, John Wiley and Sons, New York, 4th edn, 1986, Part 1, p. 13 Search PubMed.
  40. D. G. Truhlar, F. B. Brown, R. Steckler and A. D. Isaacson, The Representation and Use of Potential Energy Surfaces in the Wide Vicinity of a Reaction Path for Dynamics Calculations on Polyatomic Reactions, in The Theory of Chemical Reaction Dynamics, ed. D. C. Clary and D. Reidel, Dordrecht, Holland, 1986, p. 285 Search PubMed.
  41. S. C. Tucker and D. G. Truhlar, Dynamical Formulation of Transition State Theory: Variational Transition States and Semiclassical Tunneling, in New Theoretical Concepts for Understanding Organic Reactions, ed. J. Bertran and I. G. Csizmadia, Kluwer, Dordrecht, Netherlands, 1989, p. 291 Search PubMed.
  42. J. Zheng, J. L. Bao, R. Meana-Pañeda, S. Zhang, B. J. Lynch, J. C. Corchado, Y.-Y. Chuang, P. L. Fast, W.-P. Hu, Y.-P. Liu, G. C. Lynch, K. A. Nguyen, C. F. Jackels, A. F. Ramos, B. A. Ellingson, V. S. Melissas, J. Villà, I. Rossi, E. L. Coitino, J. Pu, T. V. Albu, R. Steckler, B. C. Garrett, A. D. Isaacson and D. G. Truhlar, Polyrate 2017 revision B, University of Minnesota, Minneapolis, 2017, http://https://comp.chem.umn.edu/polyrate/ Search PubMed.
  43. J. B. Anderson, J. Chem. Phys., 1973, 58, 4684 CrossRef CAS.
  44. J. B. Anderson, Adv. Chem. Phys., 1995, 91, 381 CrossRef CAS.
  45. A. D. Isaacson and D. G. Truhlar, J. Chem. Phys., 1982, 76, 1380 CrossRef CAS.
  46. E. B. Wilson, Jr., J. C. Decius and P. C. Cross, Molecular Vibrations, McGraw-Hill, New York, 1955 Search PubMed.
  47. P. Pechukas and E. Pollak, J. Chem. Phys., 1979, 71, 2062 CrossRef CAS.
  48. D. A. McQuarrie, Statistical Mechanics, Harper & Row, New York, 1973 Search PubMed.
  49. D. G. Truhlar and A. Kuppermann, J. Am. Chem. Soc., 1971, 93, 1840 CrossRef.
  50. B. C. Garrett, D. G. Truhlar, R. S. Grev and A. W. Magnuson, J. Phys. Chem., 1980, 84, 1730 CrossRef CAS ; erratum: 1983, 87, 4554.
  51. D. C. Chatfield, R. S. Friedman, D. G. Truhlar, B. C. Garrett and D. W. Schwenke, J. Am. Chem. Soc., 1991, 113, 486 CrossRef CAS.
  52. D. C. Chatfield, R. S. Friedman, D. G. Truhlar and D. W. Schwenke, Faraday Discuss. Chem. Soc., 1991, 91, 289 RSC.
  53. R. S. Friedman and D. G. Truhlar, Chem. Phys. Lett., 1991, 183, 539 CrossRef CAS.
  54. D. C. Chatfield, R. S. Friedman, G. C. Lynch, D. G. Truhlar and D. W. Schwenke, J. Chem. Phys., 1993, 98, 342 CrossRef CAS.
  55. R. S. Friedman, V. D. Hullinger and D. G. Truhlar, J. Phys. Chem., 1995, 99, 3184 CrossRef CAS.
  56. R. S. Friedman and D. G. Truhlar, Barrier Resonances and Chemical Reactivity, in Multiparticle Quantum Scattering with Applications to Nuclear, Atomic, and Molecular Physics, ed. D. G. Truhlar and B. Simon, Springer, New York, 1997, p. 243 Search PubMed.
  57. For some review articles, see for instances: (a) J. M. Bowman, Acc. Chem. Res., 1986, 19, 202 CrossRef CAS; (b) J. M. Bowman, S. Carter and H. Xinchuan, Int. Rev. Phys. Chem., 2003, 22, 533 CrossRef CAS; (c) J. M. Bowman, T. Carrington and H.-D. Meyer, Mol. Phys., 2008, 106, 2145 CrossRef CAS.
  58. I. M. Alecu, J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2010, 6, 2872 CrossRef CAS PubMed.
  59. J. Bloino, M. Biczysko and V. Barone, J. Chem. Theory Comput., 2012, 8, 1015 CrossRef CAS PubMed.
  60. K. M. Kuhler, D. G. Truhlar and A. D. Isaacson, J. Chem. Phys., 1996, 104, 4664 CrossRef CAS.
  61. H. H. Nielsen, Encycl. Phys., 1959, 37, 173 Search PubMed.
  62. I. M. Mills, Vibration-Rotation Structure in Asymmetric- and symmetric-Top Molecules, in Molecular Spectroscopy: Modern Research, ed. K. N. Rao and C. W. Mathews, Academic, New York, 1972, p. 115 Search PubMed.
  63. Q. Zhang, P. N. Day and D. G. Truhlar, J. Chem. Phys., 1993, 98, 4948 CrossRef CAS.
  64. D. G. Truhlar and A. D. Isaacson, J. Chem. Phys., 1991, 94, 357 CrossRef CAS.
  65. S. C. Tucker, D. G. Truhlar, B. C. Garrett and A. D. Isaacson, J. Chem. Phys., 1985, 82, 4102 CrossRef CAS.
  66. B. C. Garrett, D. G. Truhlar and G. C. Schatz, J. Am. Chem. Soc., 1986, 108, 2876 CrossRef CAS.
  67. B. C. Garrett, D. G. Truhlar, J. M. Bowman, A. F. Wagner, D. Robie, S. Arepalli, N. Presser and R. J. Gordon, J. Am. Chem. Soc., 1986, 108, 3515 CrossRef CAS.
  68. D.-h. Lu, D. Maurice and D. G. Truhlar, J. Am. Chem. Soc., 1990, 112, 6206 CrossRef CAS.
  69. D. G. Truhlar, D.-H. Lu, S. C. Tucker, X. G. Zhao, A. Gonzàlez-Lafont, T. N. Truong, D. Maurice, Y.-P. Liu and G. C. Lynch, ACS Symp. Ser., 1992, 502, 16 CrossRef CAS.
  70. A. F. Wagner and L. B. Harding, ACS Symp. Ser., 1992, 502, 48 CrossRef CAS.
  71. Y.-P. Liu, D.-h. Lu, A. Gonzàlez-Lafont, D. G. Truhlar and B. C. Garrett, J. Am. Chem. Soc., 1993, 115, 7608 Search PubMed.
  72. V. S. Melissas and D. G. Truhlar, J. Chem. Phys., 1993, 99, 3542 CrossRef CAS.
  73. T. Allison, G. C. Lynch, D. G. Truhlar and M. S. Gordon, J. Phys. Chem., 1996, 100, 13575 CrossRef CAS.
  74. H. Lin, Y. Zhao, B. A. Ellingson, J. Pu and D. G. Truhlar, J. Am. Chem. Soc., 2005, 127, 2830 CrossRef CAS PubMed.
  75. C. F. Jackels, Z. Gu and D. G. Truhlar, J. Chem. Phys., 1995, 102, 3188 CrossRef CAS.
  76. G. A. Natanson, B. C. Garrett, T. N. Truong, T. Joseph and D. G. Truhlar, J. Chem. Phys., 1991, 94, 7875 CrossRef CAS.
  77. B. C. Garrett and D. G. Truhlar, J. Phys. Chem., 1979, 83, 1079 CrossRef CAS.
  78. J. L. Bao, J. Zheng and D. G. Truhlar, J. Am. Chem. Soc., 2016, 138, 2690 CrossRef CAS PubMed.
  79. A. Fernández-Ramos, B. A. Ellingson, R. Meana-Pañeda, J. M. C. Marques and D. G. Truhlar, Theor. Chem. Acc., 2007, 118, 813 CrossRef.
  80. A. Fernández-Ramos, B. A. Ellingson, R. Meana-Pañeda, J. M. C. Marques and D. G. Truhlar, Theor. Chem. Acc., 2007, 118, 813 CrossRef.
  81. C. Eckart, Phys. Rev., 1930, 35, 1303 CrossRef CAS.
  82. B. C. Garrett and D. G. Truhlar, J. Phys. Chem., 1979, 83, 2921 CrossRef CAS.
  83. E. Wigner, Z. Phys. Chem. Abt. B, 1932, 19, 203 Search PubMed.
  84. R. T. Skodje and D. G. Truhlar, J. Phys. Chem., 1981, 85, 624 CrossRef CAS.
  85. R. T. Skodje, D. G. Truhlar and B. C. Garrett, J. Phys. Chem., 1981, 85, 3019 CrossRef CAS.
  86. Y.-P. Liu, G. C. Lynch, T. N. Truong, D.-h. Lu and D. G. Truhlar, J. Am. Chem. Soc., 1993, 115, 2408 CrossRef CAS.
  87. B. C. Garrett, D. G. Truhlar, A. F. Wagner and T. H. Dunning, J. Chem. Phys., 1983, 78, 4400 CrossRef CAS.
  88. B. C. Garrett, N. Abusalbi, D. J. Kouri and D. G. Truhlar, J. Chem. Phys., 1985, 83, 2252 CrossRef CAS.
  89. B. C. Garrett, T. Joseph, T. N. Truong and D. G. Truhlar, Chem. Phys., 1989, 136, 271 CrossRef CAS.
  90. A. Fernandez-Ramos and D. G. Truhlar, J. Chem. Phys., 2001, 114, 1491 CrossRef CAS.
  91. B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 1983, 79, 4931 CrossRef CAS.
  92. M. M. Kreevoy, D. Ostović, D. G. Truhlar and B. C. Garrett, J. Phys. Chem., 1986, 90, 3766 CrossRef CAS.
  93. R. Meana-Pañeda, D. G. Truhlar and A. Fernandez-Ramos, J. Chem. Theory Comput., 2010, 6, 3015 CrossRef PubMed.
  94. K. Gottfried and T.-M. Yan, Quantum Mechanics: Fundamentals, Springer, New York, 2nd edn, 2003, p. 216ff Search PubMed.
  95. N. T. Maitra and E. J. Heller, Phys. Rev. A: At., Mol., Opt. Phys., 1999, 61, 012107 CrossRef.
  96. J. Vaníček and E. J. Heller, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 026215 CrossRef PubMed.
  97. R. P. Bell, Trans. Faraday Soc., 1959, 55, 1 RSC.
  98. D.-h. Lu, T. N. Truong, V. S. Melissas, G. C. Lynch, Y.-P. Liu, B. C. Garrett, R. Steckler, A. D. Isaacson, S. N. Rai, G. C. Hancock, J. G. Lauderdale, T. Joseph and D. G. Truhlar, Comput. Phys. Commun., 1992, 71, 235 CrossRef CAS.
  99. D. G. Truhlar and M. S. Gordon, Science, 1990, 249, 491 CAS.
  100. R. A. Marcus, J. Chem. Phys., 1966, 45, 4493 CrossRef CAS.
  101. R. A. Marcus, J. Chem. Phys., 1968, 49, 2671 Search PubMed.
  102. R. E. Wyatt, J. Chem. Phys., 1969, 51, 3489 CrossRef CAS.
  103. D. G. Truhlar and A. Kuppermann, J. Chem. Phys., 1972, 56, 2232 CrossRef CAS.
  104. R. A. Marcus and M. E. Coltrin, J. Chem. Phys., 1977, 67, 2609 CrossRef CAS.
  105. V. K. Babamov and R. A. Marcus, J. Chem. Phys., 1978, 74, 1790 CrossRef.
  106. B. C. Garrett and D. G. Truhlar, Proc. Natl. Acad. Sci. U. S. A., 1979, 76, 4755 CrossRef CAS.
  107. W. H. Miller, N. C. Handy and J. E. Adams, J. Chem. Phys., 1980, 72, 99 CrossRef CAS.
  108. R. T. Skodje, D. G. Truhlar and B. C. Garrett, J. Chem. Phys., 1982, 77, 5955 CrossRef CAS.
  109. D. K. Bondi, J. N. L. Connor, B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 1983, 78, 5981 CrossRef CAS.
  110. D. G. Truhlar and B. C. Garrett, J. Chem. Phys., 1987, 84, 365 Search PubMed.
  111. Y. Kim, D. G. Truhlar and M. M. Kreevoy, J. Am. Chem. Soc., 1991, 113, 7837 CrossRef CAS.
  112. S. Sekušak, K. R. Liedl, B. M. Rode and A. Sabljić, J. Phys. Chem. A, 1997, 101, 4245 CrossRef.
  113. T. Loerting, K. R. Liedl and B. M. Rode, J. Chem. Phys., 1998, 109, 2672 CrossRef CAS.
  114. C. S. Tautermann, A. F. Voegele, T. Loerting and K. R. Liedl, J. Chem. Phys., 2002, 117, 1962 CrossRef CAS.
  115. C. S. Tautermann, A. F. Voegele, T. Loerting and K. R. Liedl, J. Chem. Phys., 2002, 117, 1967 CrossRef CAS.
  116. A. Fernandez-Ramos, D. G. Truhlar, J. C. Corchado and J. Espinosa-Garcia, J. Phys. Chem. A, 2002, 106, 4957 CrossRef CAS.
  117. A. Fernández-Ramos and D. G. Truhlar, J. Chem. Theory Comput., 2005, 1, 106 CrossRef PubMed.
  118. J. Pu, J. Gao and D. G. Truhlar, Chem. Rev., 2006, 106, 3140 CrossRef CAS PubMed.
  119. R. T. Skodje, D. G. Truhlar and B. C. Garrett, J. Chem. Phys., 1982, 77, 5955 CrossRef CAS.
  120. J. C. Corchado, E. L. Coitiño, Y.-Y. Chuang, P. L. Fast and D. G. Truhlar, J. Phys. Chem. A, 1998, 102, 2424 CrossRef CAS.
  121. J. Zheng and D. G. Truhlar, Faraday Discuss., 2012, 157, 59 RSC.
  122. G. A. Natanson, Chem. Phys. Lett., 1992, 190, 209 CrossRef CAS.
  123. G. A. Natanson, Mol. Phys., 1982, 46, 481 CrossRef CAS.
  124. Y.-Y. Chuang and D. G. Truhlar, J. Phys. Chem. A, 1998, 102, 242 CrossRef CAS.
  125. A. González-Lafont, J. Villà, J. M. Lluch, J. Bertrán, R. Steckler and D. G. Truhlar, J. Phys. Chem. A, 1998, 102, 3420 CrossRef.
  126. D. G. Truhlar and A. D. Isaacson, J. Chem. Phys., 1982, 77, 3516 CrossRef CAS.
  127. B. C. Garrett and D. G. Truhlar, Int. J. Quantum Chem., 1986, 29, 1463 CrossRef CAS.
  128. K. Liu, Annu. Rev. Phys. Chem., 2016, 67, 91 CrossRef CAS PubMed.
  129. D. G. Truhlar and D. A. Dixon, Direct-Mode Chemical Reactions: Classical Theories, in Atom-Molecule Collision Theory: A Guide for the Experimentalist, ed. R. B. Bernstein, Plenum Press, New York, 1979, p. 595 Search PubMed.
  130. R. Steckler, D. G. Truhlar and B. C. Garrett, J. Chem. Phys., 1986, 84, 6712 CrossRef CAS.
  131. B. C. Garrett, D. G. Truhlar, J. M. Bowman and A. F. Wagner, J. Chem. Phys., 1986, 90, 4305 CrossRef CAS.
  132. W. T. Duncan and T. N. Truong, J. Chem. Phys., 1995, 103, 9642 CrossRef CAS.
  133. T. N. Truong, J. Chem. Phys., 1995, 102, 5335 CrossRef CAS.
  134. J. Espinosa-Garcia, J. Chem. Phys., 1999, 111, 9330 CrossRef CAS.
  135. J. C. Corchado, D. G. Truhlar and J. Espinosa-Garcia, J. Chem. Phys., 2000, 112, 9375 CrossRef CAS.
  136. J. Espinosa-Garcia, J. Chem. Phys., 2000, 117, 2076 CrossRef.
  137. J. Zheng, T. Yu, E. Papajak, I. M. Alecu, S. L. Mielke and D. G. Truhlar, Phys. Chem. Chem. Phys., 2011, 13, 10885 RSC.
  138. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 1356 CrossRef CAS PubMed.
  139. B. A. Ellingson, V. A. Lynch, S. L. Mielke and D. G. Truhlar, J. Chem. Phys., 2006, 125, 084305 CrossRef PubMed.
  140. J. E. Kilpatrick and K. S. Pitzer, J. Chem. Phys., 1949, 17, 1064 CrossRef CAS.
  141. J. Zheng, S. L. Mielke, K. L. Clarkson and D. G. Truhlar, Comput. Phys. Commun., 2012, 183, 1803 CrossRef CAS.
  142. G. Voronoi, J. Reine Angew. Math., 1907, 133, 97 Search PubMed.
  143. F. Aurenhammer, ACM Comput. Surv., 1991, 23, 345 CrossRef.
  144. L. Simón-Carballido, J. L. Bao, T. V. Alves, R. Meana-Pañeda, D. G. Truhlar and A. Fernández-Ramos, J. Chem. Theory Comput., 2017, 13, 3478 CrossRef PubMed.
  145. J. L. Bao, L. Xing and D. G. Truhlar, J. Chem. Theory Comput., 2017, 13, 2511 CrossRef CAS PubMed.
  146. T. Yu, J. Zheng and D. G. Truhlar, Chem. Sci., 2011, 2, 2199 RSC.
  147. S. Winstein and N. J. Holness, J. Am. Chem. Soc., 1955, 77, 5562 CrossRef CAS.
  148. J. E. Baldwin, A. S. Raghavan, B. A. Hess and L. Smentek, J. Am. Chem. Soc., 2006, 128, 14854 CrossRef CAS PubMed.
  149. T. Yu, J. Zheng and D. G. Truhlar, J. Phys. Chem. A, 2012, 116, 297 CrossRef CAS PubMed.
  150. J. L. Bao, R. Meana-Pañeda and D. G. Truhlar, Chem. Sci., 2015, 6, 5866 RSC.
  151. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 2875 CrossRef CAS PubMed.
  152. J. L. Bao, P. Sripa and D. G. Truhlar, Phys. Chem. Chem. Phys., 2016, 18, 1032 RSC.
  153. S. N. Rai and D. G. Truhlar, J. Chem. Phys., 1983, 79, 6046 CrossRef CAS.
  154. D. M. Wardlaw and R. A. Marcus, Chem. Phys. Lett., 1984, 110, 230 CrossRef CAS.
  155. Y. Georgievskii and S. J. Klippenstein, J. Chem. Phys., 2003, 118, 5442 CrossRef CAS.
  156. S. J. Klippenstein, J. Chem. Phys., 1991, 94, 6469 CrossRef CAS.
  157. Y. Georgievskii and S. J. Klippenstein, J. Phys. Chem. A, 2003, 107, 9776 CrossRef CAS.
  158. S. J. Klippenstein, Y. Georgievskii and L. B. Harding, Phys. Chem. Chem. Phys., 2006, 8, 1133 RSC.
  159. J. Zheng, S. Zhang and D. G. Truhlar, J. Phys. Chem. A, 2008, 112, 11509 CrossRef CAS PubMed.
  160. S. E. Klippenstein, A. Wagner, S. Robertson, R. Dunbar and D. Wardlaw, Variflex – version 1.0, Argonne National laboratory, Argonne, IL, 1999 Search PubMed.
  161. J. L. Bao, X. Zhang and D. G. Truhlar, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 13606 CrossRef CAS PubMed.
  162. Y. Georgievskii and S. J. Klippenstein, J. Chem. Phys., 2005, 122, 194103 CrossRef PubMed.
  163. B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 1982, 76, 1853 CrossRef CAS.
  164. W. H. Miller, Reaction Path Hamiltonian for Polyatomic Systems: Further Developments and Applications, in Potential Energy Surfaces and Dynamics Calculations, ed. D. G. Truhlar, Plenum, New York, 1981, p. 265 Search PubMed.
  165. J. Villà, A. González-Lafont and J. M. Lluch, J. Phys. Chem., 1996, 100, 19389 CrossRef.
  166. L. Masgrau, A. González-Lafont and J. M. Lluch, J. Comput. Chem., 1999, 20, 1685 CrossRef CAS.
  167. A. Fernández-Ramos and A. J. C. Varandas, J. Phys. Chem. A, 2002, 106, 4077 CrossRef.
  168. E. Martínez-Núñez, A. Fernández-Ramos, S. A. Vázquez and M. A. Ríos, Chem. Phys. Lett., 2002, 360, 50 CrossRef.
  169. L. Masgrau, A. González-Lafont and J. M. Lluch, J. Phys. Chem. A, 2002, 106, 11760 CrossRef CAS.
  170. L. Masgrau, A. González-Lafont and J. M. Lluch, J. Phys. Chem. A, 1999, 103, 1044 CrossRef CAS.
  171. L. Masgrau, A. González-Lafont and J. M. Lluch, J. Chem. Phys., 2001, 114, 2154 CrossRef CAS.
  172. J. M. Ramírez-Anguita, A. González-Lafont and J. M. Lluch, J. Comput. Chem., 2011, 32, 2104 CrossRef PubMed.
  173. L. Masgrau, A. González-Lafont and J. M. Lluch, J. Phys. Chem. A, 2003, 107, 4490 CrossRef CAS.
  174. H. Ma, X. Liu, W. Bian, L. Meng and S. Zheng, ChemPhysChem, 2006, 7, 1786 CrossRef CAS PubMed.
  175. A. González-Lafont and J. M. Lluch, THEOCHEM, 2004, 709, 35 CrossRef.
  176. J. M. Ramírez-Anguita, R. Gelabert, A. González-Lafont, M. Moreno and J. M. Lluch, Theor. Chem. Acc., 2011, 128, 569 CrossRef.
  177. Y. Peng, Z. Jiang and J. Chen, J. Phys. Chem. A, 2017, 121, 2209 CrossRef CAS PubMed.
  178. (a) E. W. Montroll and K. E. Shuler, Adv. Chem. Phys., 1958, 1, 361 CrossRef CAS; (b) R. Zwanzig, Physica, 1964, 30, 1109 CrossRef; (c) I. Oppenheim and K. E. Shuler, Phys. Rev., 1965, 138, B1007 CrossRef; (d) N. S. Snider, J. Chem. Phys., 1965, 42, 548 CrossRef CAS; (e) R. K. Boyd, J. Chem. Phys., 1974, 60, 1214 CrossRef CAS; (f) H. O. Pritchard, The Quantum Theory of Unimolecular Reactions, Cambridge University Press, Cambridge, 1984 Search PubMed; (g) J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2006, 110, 10528 CrossRef CAS PubMed.
  179. R. G. Gilbert and S. C. Smith, Theory of Unimolecular and Recombination Reactions, Blackwell Science, Oxford, 1990 Search PubMed.
  180. S. J. Klippenstein, V. Pande and D. G. Truhlar, J. Am. Chem. Soc., 2014, 136, 528 CrossRef CAS PubMed.
  181. K. A. Holbrook, M. J. Pilling and S. H. Robertson, Unimolecular Reactions, John Wiley & Sons, New York, 2nd edn, 1996 Search PubMed.
  182. J. R. Barker and D. M. Golden, Chem. Rev., 2003, 103, 4577 CrossRef CAS PubMed.
  183. J. W. Allen, C. F. Goldsmith and W. H. Green, Phys. Chem. Chem. Phys., 2012, 14, 1131 RSC.
  184. D. R. Glowacki, C.-H. Liang, C. Morley, M. J. Pilling and S. H. Robertson, J. Phys. Chem. A, 2012, 116, 9545 CrossRef CAS PubMed.
  185. S. H. Robertson, D. R. Glowacki, C.-H. Liang, C. Morley and M. J. Pilling, MESMER (Master Equation Solver for Multi-Energy Well Reactions), 2008, http://sourceforge.net/projects/mesmer Search PubMed.
  186. J. R. Barker, Int. J. Chem. Kinet., 2001, 33, 232 CrossRef CAS.
  187. J. W. Allen, R. W. Ashcraft, G. J. Beran, C. F. Goldsmith, M. R. Harper, A. Jalan, G. R. Magoon, D. M. Matheu, S. S. Merchant, J. D. Mo, S. Petway, S. Ruman, S. Sharma, K. M. Van Geem, J. Song, J. Wen, R. H. West, A. Wong, H.-W. Wong, P. E. Yelvington, J. Yu and W. H. Green, RMG (Reaction Mechanism Generator) version 3.2, 2010, http://rmg. sourceforge.net/ Search PubMed.
  188. For examples of some earlier work, see: (a) B. S. Rabinovitch, D. H. Dills, W. H. McLain and J. H. Current, J. Chem. Phys., 1960, 32, 493 CrossRef CAS; (b) R. E. Harrington, B. S. Rabinovitch and R. W. Diesen, J. Chem. Phys., 1960, 32, 1245 CrossRef CAS; (c) R. E. Harrington, B. S. Rabinovitch and M. R. Hoare, J. Chem. Phys., 1960, 33, 744 CrossRef CAS.
  189. A. M. Dean, J. Phys. Chem., 1985, 89, 4600 CrossRef CAS.
  190. J. R. Barker, J. Chem. Phys., 1980, 72, 3686 CrossRef CAS.
  191. J. R. Barker, Chem. Phys., 1983, 77, 301 CrossRef CAS.
  192. P. R. Westmoreland, J. B. Howard, J. P. Longwell and A. M. Dean, AIChE J., 1986, 32, 1971 CrossRef CAS.
  193. J. W. Bozzelli and A. M. Dean, J. Phys. Chem., 1990, 94, 3313 CrossRef CAS.
  194. J. W. Bozzelli, A. Y. Chang and A. M. Dean, Int. J. Chem. Kinet., 1997, 29, 161 CrossRef CAS.
  195. A. Y. Chang, J. W. Bozzelli and A. M. Dean, Z. Phys. Chem., 2000, 214, 1533 CrossRef CAS.
  196. H. H. Carstensen and A. M. Dean, Compr. Chem. Kinet., 2007, 42, 101 Search PubMed.
  197. C. Y. Sheng, J. W. Bozzelli, A. M. Dean and A. Y. Chang, J. Phys. Chem. A, 2002, 106, 7276 CrossRef CAS.
  198. J. L. Bao, X. Zhang and D. G. Truhlar, Phys. Chem. Chem. Phys., 2016, 18, 16659 RSC.
  199. J. L. Bao and D. G. Truhlar, Phys. Chem. Chem. Phys., 2016, 18, 10097 RSC.
  200. L. S. Kassel, J. Phys. Chem., 1928, 32, 1065 CrossRef CAS.
  201. L. S. Kassel, Proc. Natl. Acad. Sci. U. S. A., 1928, 14, 23 CrossRef CAS.
  202. J. Zheng and D. G. Truhlar, Phys. Chem. Chem. Phys., 2010, 12, 7782 RSC.
  203. D. G. Truhlar, J. Chem. Educ., 1978, 55, 309 CrossRef CAS.
  204. A. Dollet, S. de Persis and F. Teyssandier, Phys. Chem. Chem. Phys., 2004, 6, 1203 RSC.
  205. J. R. Barker and N. F. Ortiz, Int. J. Chem. Kinet., 2001, 33, 246 CrossRef CAS.
  206. J. Troe, J. Chem. Phys., 1977, 66, 4745 CrossRef CAS.
  207. J. Troe, J. Chem. Phys., 1977, 66, 4758 CrossRef CAS.
  208. J. Troe, Ber. Bunsen-Ges., 1983, 87, 161 CrossRef CAS.
  209. R. G. Gilbert, K. Luther and J. Troe, Ber. Bunsen-Ges., 1983, 87, 169 CrossRef CAS.
  210. H. Hippler, H. J. Wendelken and J. Troe, J. Chem. Phys., 1983, 78, 6709 CrossRef CAS.
  211. D. C. Tardy and B. S. Rabinovitch, J. Phys. Chem., 1985, 89, 2442 CrossRef CAS.
  212. J. Troe, J. Phys. Chem., 1979, 83, 114 CrossRef CAS.
  213. J. Troe and V. G. Ushakov, J. Chem. Phys., 2011, 135, 054304 CrossRef CAS PubMed.
  214. A. J. Stace and J. N. Murrell, J. Chem. Phys., 1978, 68, 3028 CrossRef CAS.
  215. A. W. Jasper, J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2013, 117, 12243 CrossRef CAS PubMed.
  216. A. W. Jasper, K. M. Pelzer, J. A. Miller, E. Kamarchik, L. B. Harding and S. J. Klippenstein, Science, 2014, 346, 1212 CrossRef CAS PubMed.
  217. A. W. Jasper, C. M. Oana and J. A. Miller, Proc. Combust. Inst., 2015, 35, 197 CrossRef CAS.
  218. S. J. Klippenstein, Proc. Combust. Inst., 2017, 36, 77 CrossRef CAS.
  219. D. Rapp and T. Kassal, Chem. Rev., 1969, 69, 61 CrossRef CAS.
  220. (a) S. H. Luu and J. Troe, Ber. Bunsen-Ges. Phys. Chem., 1974, 78, 766 CAS; (b) A. A. Shestov, S. A. Kostina and V. D. Knyazev, Proc. Combust. Inst., 2005, 30, 975 CrossRef.
  221. A. W. Jasper and J. A. Miller, J. Phys. Chem. A, 2009, 113, 5612 CrossRef CAS PubMed.
  222. J. R. Barker, J. Phys. Chem., 1984, 88, 11 CrossRef CAS.
  223. M. L. Yerram, J. D. Brenner, K. D. King and J. R. Barker, J. Phys. Chem., 1990, 94, 6341 CrossRef CAS.
  224. I. Oref and D. C. Tardy, Chem. Rev., 1990, 90, 1407 CrossRef CAS.
  225. G. W. Flynn, C. S. Parmenter and A. M. Wodtke, J. Phys. Chem., 1996, 100, 12817 CrossRef CAS.
  226. T. Lenzer, K. Luther, K. Reihs and A. C. Symonds, J. Chem. Phys., 2000, 112, 4090 CrossRef CAS.
  227. J. R. Barker, L. M. Yoder and K. D. King, J. Phys. Chem. A, 2001, 105, 796 CrossRef CAS.
  228. J. Du, L. Yuan, S. Hsieh, F. Lin and A. S. Mullin, J. Phys. Chem. A, 2008, 112, 9396 CrossRef CAS PubMed.
  229. D. C. Tardy and B. S. Rabinovitch, Chem. Rev., 1977, 77, 369 CrossRef CAS.
  230. J. Troe, Ber. Bunsen-Ges. Phys. Chem., 1983, 87, 161 CrossRef CAS.
  231. R. G. Gilbert, K. Luther and J. Troe, Ber. Bunsen-Ges. Phys. Chem., 1983, 87, 169 CrossRef CAS.
  232. B. Long, J. L. Bao and D. G. Truhlar, Phys. Chem. Chem. Phys., 2017, 19, 8091 RSC.
  233. W. H. Miller, J. Chem. Phys., 1974, 61, 1823 CrossRef CAS.
  234. P. Pechukas, Ber. Bunsen-Ges. Phys. Chem., 1982, 86, 372 CrossRef CAS.
  235. R. P. Bell, W. H. Sachs and R. L. Tranter, Trans. Faraday Soc., 1971, 67, 1995 RSC.
  236. J. Pang, S. Hay, N. S. Scrutton and M. J. Sutcliffe, J. Am. Chem. Soc., 2008, 130, 7092 CrossRef CAS PubMed.
  237. A. Wu and J. M. Mayer, J. Am. Chem. Soc., 2008, 130, 14745 CrossRef CAS PubMed.
  238. H. S. Johnston, Gas Phase Reaction Rate Theory, Ronald Press, New York, 1966, pp. 133–135 Search PubMed.
  239. B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 1980, 72, 3460 CrossRef CAS.
  240. J. Pu and D. G. Truhlar, J. Chem. Phys., 2002, 117, 1479 CrossRef CAS.
  241. R. A. Marcus, Faraday Discuss. Chem. Soc., 1973, 55, 9 RSC.
  242. H. B. Wang, G. H. Peslherbe and W. L. Hase, J. Am. Chem. Soc., 1994, 116, 9644 CrossRef CAS.
  243. H. B. Wang, L. Zhu and W. L. Hase, J. Phys. Chem., 1994, 98, 1608 CrossRef CAS.
  244. J. Xie and W. L. Hase, Science, 2016, 352, 32 CrossRef CAS PubMed.
  245. J. O. Hirschfelder and E. Wigner, J. Chem. Phys., 1939, 7, 616 CrossRef CAS.
  246. W. H. Miller, J. Chem. Phys., 1976, 65, 2216 CrossRef CAS.
  247. W.-P. Hu and D. G. Truhlar, J. Am. Chem. Soc., 1995, 117, 10726 CrossRef CAS.
  248. W.-P. Hu and D. G. Truhlar, J. Am. Chem. Soc., 1996, 118, 860 CrossRef CAS.
  249. J. A. Miller, J. Chem. Phys., 1981, 74, 5120 CrossRef CAS.
  250. W. L. Hase, S. L. Mondro, R. J. Duchovic and D. M. Hirst, J. Am. Chem. Soc., 1987, 109, 2916 CrossRef CAS.
  251. R. A. Marcus, Faraday Discuss. Chem. Soc., 1973, 55, 379 Search PubMed.
  252. B. R. Beno, K. N. Houk and D. A. Singleton, J. Am. Chem. Soc., 1996, 118, 9984 CrossRef CAS.
  253. Y. Nieves-Quinones and D. A. Singleton, J. Am. Chem. Soc., 2016, 138, 15167 CrossRef CAS PubMed.
  254. H. Kurouchi, I. L. andujar-De Sanctis and D. A. Singleton, J. Am. Chem. Soc., 2016, 138, 14534 CrossRef CAS PubMed.
  255. C. Doubleday, M. Boguslav, C. Howell, S. D. Korotkin and D. Shaked, J. Am. Chem. Soc., 2016, 138, 7476 CrossRef CAS PubMed.
  256. J. Rehbein and B. K. Carpenter, Phys. Chem. Chem. Phys., 2011, 13, 20906 RSC.
  257. J. Zheng, E. Papajak and D. G. Truhlar, J. Am. Chem. Soc., 2009, 131, 15754 CrossRef CAS PubMed.
  258. R. V. Serauskas and E. W. Schlag, J. Chem. Phys., 1966, 45, 3706 CrossRef CAS.
  259. D. G. Truhlar, J. Am. Chem. Soc., 1975, 97, 6310 CrossRef CAS.
  260. P. Pechukas and J. C. Light, J. Chem. Phys., 1965, 42, 3281 CrossRef CAS.
  261. E. E. Nikitin, Theor. Exp. Chem., 1965, 1, 275 CrossRef.
  262. D. R. Glowacki, C. H. Liang, S. P. Marsden, J. N. Harvey and M. J. Pilling, J. Am. Chem. Soc., 2010, 132, 13621 CrossRef CAS PubMed.
  263. D. G. Truhlar, Variational Transition-State Theory and Multidimensional Tunneling for Simple and Complex Reactions in the Gas Phase, Solids, Liquids and Enzymes, in Isotope Effects in Chemistry and Biology, ed. A. Kohen and H.-H. Limbach, Marcel Dekker, Inc., New York, 2006, p. 579 Search PubMed.
  264. Y. Kim, J. R. Mohrig and D. G. Truhlar, J. Am. Chem. Soc., 2010, 132, 11071 CrossRef CAS PubMed.
  265. C. J. Cramer and D. G. Truhlar, The Thermodynamics of Solvation and the Treatment of Equilibrium and Nonequilibrium Solvation Effects by Models Based on Collective Solvent Coordinates, in Free Energy Calculations in Rational Drug Design, ed. M. R. Reddy and M. D. Erion, Kluwer Academic/Plenum, New York, 2001, p. 63 Search PubMed.
  266. D. G. Truhlar, Molecular-Scale Modeling of Reactions and Solvation, in First International Conference on Foundations of Molecular Modeling and Simulation, ed. P. T. Cummings, P. R. Westmoreland and B. Carnahan, American Institute of Chemical Engineers, New York, Symposium Series, 2001, vol. 97, p. 71 Search PubMed.
  267. C. J. Cramer and D. G. Truhlar, Development and Biological Applications of Quantum Mechanical Continuum Solvation Models in Quantitative Treatments of Solute/Solvent Interactions, ed. P. Politzer and J. S. Murray, Elsevier, Amsterdam, 1994, p. 83 Search PubMed.
  268. J.-L. Rivail and D. Rinaldi, Liquid-State Quantum Chemistry: Computational Applications of the Polarizable Continuum Models in Computational Chemistry: Reviews of Current Trends, ed. J. Leszczynynski, World Scientific, Singapore, 1996, p. 139 Search PubMed.
  269. C. J. Cramer and D. G. Truhlar, Chem. Rev., 1999, 99, 2161 CrossRef CAS PubMed.
  270. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999 CrossRef CAS PubMed.
  271. C. J. Cramer and D. G. Truhlar, SMx Continuum Models for Condensed Phases in Trends and Perspectives in Modern Computational Science, Lecture Series on Computer and Computational Sciences, ed. G. Maroulis and T. E. Simos, Brill/VSP, Leiden, 2006, vol. 6, p. 112 Search PubMed.
  272. C. J. Cramer and D. G. Truhlar, Acc. Chem. Res., 2008, 41, 760 CrossRef CAS PubMed.
  273. D. G. Truhlar and J. R. Pliego Jr., Transition State Theory and Chemical Reaction Dynamics in Solution in Continuum Solvation Models in Chemical Physics: From Theory to Applications, ed. B. Mennucci and R. Cammi, Wiley, Chichester, 2008, p. 338 Search PubMed.
  274. C. P. Kelly, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. A, 2006, 110, 2493 CrossRef CAS PubMed.
  275. C. P. Kelly, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2006, 110, 16066 CrossRef CAS PubMed.
  276. A. V. Marenich, W. Ding, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. Lett., 2012, 3, 1437 CrossRef CAS PubMed.
  277. A. Ben-Naim, Statistical Thermodynamics for Chemists and Biochemists, Plenum, New York, 1992 Search PubMed.
  278. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378 CrossRef CAS PubMed.
  279. R. F. Ribeiro, A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2011, 115, 14556 CrossRef CAS PubMed.
  280. D. G. Truhlar, Y.-P. Liu, G. K. Schenter and B. C. Garrett, J. Phys. Chem., 1994, 98, 8396 CrossRef CAS.
  281. Y.-Y. Chuang, C. J. Cramer and D. G. Truhlar, Int. J. Quantum Chem., 1998, 70, 887 CrossRef CAS.
  282. G. K. Schenter, B. C. Garrett and D. G. Truhlar, J. Chem. Phys., 2003, 119, 5828 CrossRef CAS.
  283. G. van der Zwan and J. T. Hynes, J. Chem. Phys., 1983, 78, 4174 CrossRef CAS.
  284. H. A. Kramers, Physica, 1940, 7, 284 CrossRef CAS.
  285. R. F. Grote and J. T. Hynes, J. Chem. Phys., 1980, 73, 2715 CrossRef CAS.
  286. S. C. Tucker, Variational Transition State theory in Condensed Phases, in New Trends in Kramers’ Reaction Rate Theory, ed. P. Talkner and P. Hänggi, Kluwer Academic, Dordrecht, 1995, p. 5 Search PubMed.
  287. D. G. Truhlar and B. C. Garrett, J. Phys. Chem. B, 2000, 104, 1069 CrossRef CAS.
  288. N. G. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amterdam, 1981 Search PubMed.
  289. B. C. Garrett and G. K. Schenter, Int. Rev. Phys. Chem., 1994, 13, 263 CrossRef CAS.
  290. B. C. Garrett and G. K. Schenter, ACS Symp. Ser., 1994, 568, 122 CrossRef CAS.
  291. Y.-Y. Chuang and D. G. Truhlar, J. Am. Chem. Soc., 1999, 121, 10157 CrossRef CAS.
  292. A. F. Voter and J. D. Doll, J. Chem. Phys., 1984, 80, 5832 CrossRef CAS.
  293. J. D. Doll, J. Chem. Phys., 1980, 73, 2760 CrossRef CAS.
  294. J. D. Doll, J. Chem. Phys., 1981, 74, 1074 CrossRef CAS.
  295. N. Metropolis, A. Rosenbluth, M. N. Rosenbluth, A. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087 CrossRef CAS.
  296. J. G. Lauderdale and D. G. Truhlar, Surf. Sci., 1985, 164, 558 CrossRef CAS.
  297. J. G. Lauderdale and D. G. Truhlar, J. Chem. Phys., 1986, 84, 1843 CrossRef CAS.
  298. S. E. Wonchoba and D. G. Truhlar, J. Chem. Phys., 1993, 99, 9637 CrossRef CAS.
  299. S. E. Wonchoba, W.-P. Hu and D. G. Truhlar, Reaction Path Approach to Dynamics at a Gas-Solid Interface: Quantum Tunneling Effects for an Adatom on a Non-Rigid Metallic Surface, in Theoretical and Computational Approaches to Interface Phenomena, ed. H. L. Sellers and J. T. Golab, Plenum, New York, 1994, p. 1 Search PubMed.
  300. S. E. Wonchoba, W.-P. Hu and D. G. Truhlar, Phys. Rev. B: Condens. Matter Mater. Phys., 1995, 51, 9985 CrossRef CAS.
  301. J. G. Lauderdale and D. G. Truhlar, J. Am. Chem. Soc., 1985, 107, 4590 CrossRef CAS.
  302. T. N. Truong, D. G. Truhlar, J. R. Chelikowsky and M. Y. Chou, J. Phys. Chem., 1990, 94, 1973 CrossRef CAS.
  303. T. N. Truong and D. G. Truhlar, J. Chem. Phys., 1990, 93, 2125 CrossRef CAS.
  304. T. N. Truong and D. G. Truhlar, J. Phys. Chem., 1990, 94, 8262 CrossRef CAS.
  305. T.-S. Lin and R. Gomer, Surf. Sci., 1991, 225, 41 CrossRef.
  306. X. D. Zhu, A. Lee, A. Wong and U. Linke, Phys. Rev. Lett., 1992, 68, 1862 CrossRef CAS PubMed.
  307. A. Lee, X. D. Zhu, L. Deng and U. Linke, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 15472 CrossRef CAS.
  308. P. S. Zuev, R. S. Sheridan, T. V. Albu, D. G. Truhlar, D. A. Hrovat and W. T. Borden, Science, 2003, 299, 867 CrossRef CAS PubMed.
  309. C. Alhambra, J. Corchado, M. L. Sánchez, M. Garcia-Viloca, J. Gao and D. G. Truhlar, J. Phys. Chem. B, 2001, 105, 11326 CrossRef CAS.
  310. L. Masgrau and D. G. Truhlar, Acc. Chem. Res., 2015, 48, 431 CrossRef CAS PubMed.
  311. D. G. Truhlar, J. Gao, C. Alhambra, M. Garcia-Viloca, J. Corchado, M. L. Sanchez and J. Villa, Acc. Chem. Res., 2002, 35, 341 CrossRef CAS PubMed.
  312. M. Garcia-Viloca, C. Alhambra, D. G. Truhlar and J. Gao, J. Chem. Phys., 2001, 114, 9953 CrossRef CAS.
  313. M. Garcia-Viloca, D. G. Truhlar and J. Gao, Biochemistry, 2003, 42, 13558 CrossRef CAS PubMed.
  314. D. G. Truhlar, J. Phys. Org. Chem., 2010, 23, 660 CrossRef CAS.
  315. C. A. Parr and D. G. Truhlar, J. Phys. Chem., 1971, 75, 1844 CrossRef.
  316. M. Monge-Palacios, C. Rangel, J. C. Corchado and J. Espinosa-García, Int. J. Quantum Chem., 2011, 112, 1887 CrossRef.
  317. D. G. Truhlar, R. Steckler and M. S. Gordon, Chem. Rev., 1987, 87, 217 CrossRef CAS.
  318. T. Albu, J. C. Corchado and D. G. Truhlar, J. Phys. Chem. A, 2001, 105, 8465 CrossRef CAS.
  319. T. V. Albu, J. Espinosa-García and D. G. Truhlar, Chem. Rev., 2007, 107, 5101 CrossRef CAS PubMed.
  320. E. Gonzalez-Lavado, J. C. Corchado, Y. V. Suleimanov, W. H. Green and J. Espinosa-Garcia, J. Phys. Chem. A, 2014, 118, 3243 CrossRef CAS PubMed.
  321. M. Monge-Palacios, C. Rangel and J. Espinosa-Garcia, J. Chem. Phys., 2013, 138, 084305 CrossRef CAS PubMed.
  322. M. Monge-Palacios, C. Rangel and J. Espinosa-Garcia, Phys. Chem. Chem. Phys., 2016, 18, 16941 RSC.
  323. J. Espinosa-García, J. L. Bravo and C. Rangel, J. Phys. Chem. A, 2007, 111, 2761 CrossRef PubMed.
  324. J. C. Corchado and J. Espinosa-García, J. Chem. Phys., 1996, 105, 3152 CrossRef CAS.
  325. K. K. Baldridge, M. S. Gordon, R. Steckler and D. G. Truhlar, J. Phys. Chem., 1989, 93, 5107 CrossRef CAS.
  326. D. G. Truhlar, Direct Dynamics Method for the Calculation of Reaction Rates, in The Reaction Path in Chemistry: Current Approaches and Perspectives, ed. D. Heidrich, Kluwer, Dordrecht, 1995, p. 229 Search PubMed.
  327. J. Zheng, P. Seal and D. G. Truhlar, Chem. Sci., 2013, 4, 200 RSC.
  328. J. Zheng, R. Meana-Pañeda and D. G. Truhlar, J. Am. Chem. Soc., 2014, 136, 5150 CrossRef CAS PubMed.
  329. J. A. Sansón, M.-L. Sánchez and J. C. Corchado, J. Phys. Chem. A, 2006, 110, 589 CrossRef PubMed.
  330. A. D. Isaacson, J. Phys. Chem. A, 2006, 110, 379 CrossRef CAS PubMed.
  331. J. Park, Z. F. Xu and M. C. Lin, J. Chem. Phys., 2003, 118, 9990 CrossRef CAS.
  332. T. Joseph, R. Steckler and D. G. Truhlar, J. Chem. Phys., 1987, 87, 7036 CrossRef CAS.
  333. R. Meana-Pañeda and A. Fernández-Ramos, J. Chem. Phys., 2014, 140, 174303 CrossRef PubMed.
  334. P. Seal, E. Papajak and D. G. Truhlar, J. Phys. Chem. Lett., 2012, 3, 264 CrossRef CAS PubMed.
  335. C.-W. Zhou, J. M. Simmie and H. J. Curran, Int. J. Chem. Kinet., 2012, 44, 155 CrossRef CAS.
  336. I. M. Alecu and D. G. Truhlar, J. Phys. Chem. A, 2011, 115, 14599 CrossRef CAS PubMed.
  337. I. M. Alecu, J. Zheng, E. Papajak, T. Yu and D. G. Truhlar, J. Phys. Chem. A, 2012, 116, 12206 CrossRef CAS PubMed.
  338. P. Seal, G. Oyedepo and D. G. Truhlar, J. Phys. Chem. A, 2013, 117, 275 CrossRef CAS PubMed.
  339. J. Zheng, G. Oyedepo and D. G. Truhlar, J. Phys. Chem. A, 2015, 119, 12182 CrossRef CAS PubMed.
  340. A. D. Isaacson, M. T. Sund, S. N. Rai and D. G. Truhlar, J. Chem. Phys., 1985, 82, 1338 CrossRef CAS.
  341. M. Akbar Ali and J. R. Barker, J. Phys. Chem. A, 2015, 119, 7578 CrossRef CAS PubMed.
  342. X. Jia, Y. Liu, J. Sun, H. Sun, Z. Su, X. Pan and R. Wang, J. Phys. Chem. A, 2010, 114, 417 CrossRef CAS PubMed.
  343. S. R. Sellevåg, G. Nyman and C. J. Nielsen, J. Phys. Chem. A, 2006, 110, 141 CrossRef PubMed.
  344. L. Yang, J.-Y. Liu, S.-Q. Wan and Z.-S. Li, J. Comput. Chem., 2009, 30, 565 CrossRef CAS PubMed.
  345. L. Sandhiya, P. Kolandaivel and K. Senthilkumar, J. Phys. Chem. A, 2013, 117, 4611 CrossRef CAS PubMed.
  346. M. R. Dash and B. Rajakumar, Chem. Phys. Lett., 2014, 597, 75 CrossRef CAS.
  347. S. Sekušak, K. R. Liedl, B. M. Rode and A. Sabljić, J. Phys. Chem. A, 1997, 101, 4245 CrossRef.
  348. L. Yang, J.-Y. Liu and Z.-S. Li, J. Phys. Chem. A, 2008, 112, 6364 CrossRef CAS PubMed.
  349. F.-Y. Bai, G. Sun, X. Wang, Y.-Q. Sun, R.-S. Wang and X.-M. Pan, J. Phys. Chem. A, 2015, 119, 1256 CrossRef CAS PubMed.
  350. G. Srinivasulu and B. Rajakumar, J. Phys. Chem. A, 2013, 117, 4534 CrossRef CAS PubMed.
  351. L. Yang, J.-Y. Liu, L. Wang, H.-Q. He, Y. Wang and Z.-S. Li, J. Phys. Chem. A, 2008, 29, 550 CAS.
  352. S. Sekušak, M. G. Cory, R. J. Bartlett and A. Sabljić, J. Phys. Chem. A, 1999, 103, 11394 CrossRef.
  353. X. Xu, T. Yu, E. Papajak and D. G. Truhlar, J. Phys. Chem. A, 2012, 116, 10480 CrossRef CAS PubMed.
  354. C. Zavala-Oseguera, J. R. Alvarez-Idaboy, G. Merine and A. Galano, J. Phys. Chem. A, 2009, 113, 13913 CrossRef CAS PubMed.
  355. W.-P. Hu, I. Rossi, J. C. Corchado and D. G. Truhlar, J. Phys. Chem. A, 1997, 101, 6911 CrossRef CAS.
  356. D. G. Truhlar, K. Runge and B. C. Garrett, Variational Transition State Theory and Tunneling Calculations of Potential Energy Surface Effects on the Reaction of O(3P) with H2, in Twentieth Symposium (International) on Combustion, Combustion Institute, Pittsburgh, 1984, p. 585 Search PubMed.
  357. G. Da Silva and J. W. Bozzelli, Phys. Chem. Chem. Phys., 2012, 14, 16143 RSC.
  358. R. F. K. Spada, L. F. A. Ferrão, R. J. Rocha, K. Iha, J. A. F. F. Rocco, O. Roberto-Neto, H. Lischka and F. B. C. Machado, J. Phys. Chem. A, 2015, 119, 1628 CrossRef CAS PubMed.
  359. B. C. Garrett and D. G. Truhlar, Int. J. Quantum Chem., 1987, 31, 17 CrossRef CAS.
  360. A. G. S. De Oliveira-Filho, F. R. Ornellas, K. A. Peterson and S. L. Mielke, J. Phys. Chem. A, 2013, 117, 12703 CrossRef CAS PubMed.
  361. J. C. Corchado, J. Espinosa-Garcia, O. Roberto-Neto, Y.-Y. Chuang and D. G. Truhlar, J. Phys. Chem. A, 1998, 102, 4899 CrossRef CAS.
  362. R. Meana-Pañeda, X. Xu, H. Ma and D. G. Truhlar, J. Phys. Chem. A, 2017, 121, 1693 CrossRef PubMed.
  363. A. Jalan, I. M. Alecu, R. Meana-Pañeda, J. Aguilera-Iparraguirre, K. R. Yang, S. S. Merchant, D. G. Truhlar and W. H. Green, J. Am. Chem. Soc., 2013, 135, 11100 CrossRef CAS PubMed.
  364. B. C. Garrett and D. G. Truhlar, J. Am. Chem. Soc., 1979, 101, 5207 CrossRef CAS.
  365. T. V. Mai, A. Ratkiewicz, M. V. Duong and L. K. Huynh, Chem. Phys. Lett., 2016, 646, 102 CrossRef CAS.
  366. M. L. Coote, J. Phys. Chem. A, 2005, 109, 1230 CrossRef CAS PubMed.
  367. N. Kungwan and T. N. Truong, J. Phys. Chem. A, 2005, 109, 7742 CrossRef CAS PubMed.
  368. A. Datta, D. A. Hrovat and W. T. Borden, J. Am. Chem. Soc., 2008, 130, 6684 CrossRef CAS PubMed.
  369. S. Sun, S. Cheng and H. Zhang, Theor. Chem. Acc., 2016, 135, 154 CrossRef.
  370. Y. Ge, M. S. Gordon, F. Battaglia and R. O. Fox, J. Phys. Chem. A, 2010, 114, 2384 CrossRef CAS PubMed.
  371. J. Zheng and D. G. Truhlar, J. Phys. Chem. A, 2009, 113, 11919 CrossRef CAS PubMed.
  372. J. Zheng and D. G. Truhlar, Phys. Chem. Chem. Phys., 2010, 12, 7782 RSC.
  373. K. T. Kuwata, M. R. Hermes, M. J. Carlson and C. K. Zogg, J. Phys. Chem. A, 2010, 114, 9192 CrossRef CAS PubMed.
  374. X. Xu, E. Papajak, J. Zheng and D. G. Truhlar, Phys. Chem. Chem. Phys., 2012, 14, 4204 RSC.
  375. A. Ratkiewicz, B. Bankiewicz and T. N. Truong, Phys. Chem. Chem. Phys., 2010, 12, 10988 RSC.
  376. K. T. Kuwata, T. S. Dibble, E. Sliz and E. B. Petersen, J. Phys. Chem. A, 2007, 111, 5032 CrossRef CAS PubMed.
  377. F. Zhang and T. S. Dibble, Phys. Chem. Chem. Phys., 2011, 13, 17969 RSC.
  378. Y. Guan, Y. Li, L. Zhao, Y. Song, H. Ma and J. Song, Comput. Theor. Chem., 2016, 1089, 43 CrossRef CAS.
  379. Y. Guan, J. Gao, Y. Song, Y. Li, H. Ma and J. Song, J. Phys. Chem. A, 2017, 121, 1121 CrossRef CAS PubMed.
  380. Y. Sakai, H. Ando, H. K. Chakravarty, H. Pitsch and R. X. Fernandes, Proc. Combust. Inst., 2015, 35, 161 CrossRef CAS.
  381. J. A. Miller and S. J. Klippenstein, Phys. Chem. Chem. Phys., 2004, 6, 1192 RSC.
  382. A. W. Jasper, S. J. Klippenstein, L. B. Harding and B. Ruscic, J. Phys. Chem. A, 2007, 111, 3932 CrossRef CAS PubMed.
  383. S. J. Klippenstein, J. A. Miller and A. W. Jasper, J. Phys. Chem. A, 2015, 119, 7780 CrossRef CAS PubMed.
  384. E. E. Greenwald, S. W. North, Y. Georgievskii and S. J. Klippenstein, J. Phys. Chem. A, 2005, 109, 6031 CrossRef CAS PubMed.
  385. B. Long, J. L. Bao and D. G. Truhlar, J. Am. Chem. Soc., 2016, 138, 14409 CrossRef CAS PubMed.
  386. (a) M. A. Blitz, R. J. Salter, D. E. Heard and P. W. Seakins, J. Phys. Chem. A, 2017, 121, 3175 CrossRef CAS PubMed; (b) M. A. Blitz, R. J. Salter, D. E. Heard and P. W. Seakins, J. Phys. Chem. A, 2017, 121, 3184 CrossRef CAS PubMed.
  387. A. Cartoni, D. Catone, P. Bolognesi, M. Satta, P. Markus and L. Avaldi, Chem. – Eur. J., 2017, 23, 6772 CrossRef CAS PubMed.
  388. B. Long, W.-J. Zhang, X.-F. Tan, Z.-W. Long, Y.-B. Wang and D.-S. Ren, J. Phys. Chem. A, 2011, 115, 1350 CrossRef CAS PubMed.
  389. B. Long, X.-F. Tan, Y.-B. Wang, J. Li, D.-S. Ren and W.-J. Zhang, ChemistrySelect, 2016, 1, 1421 CrossRef CAS.
  390. B. A. Ellingson and D. G. Truhlar, J. Am. Chem. Soc., 2007, 129, 12765 CrossRef CAS PubMed.
  391. T. Loerting, K. R. Liedl and B. M. Rode, J. Chem. Phys., 1998, 109, 2672 CrossRef CAS.
  392. For instances, (a) B. Long, X.-F. Tan, J. L. Bao, D.-M. Wang and Z.-W. Long, Int. J. Chem. Kinet., 2016, 49, 130 CrossRef; (b) B. Long, X.-F. Tan, C.-R. Chang, W.-X. Zhao, Z.-W. Long, D.-S. Ren and W.-J. Zhang, J. Phys. Chem. A, 2013, 117, 5106 CrossRef CAS PubMed; (c) B. Long, C.-R. Chang, Z.-W. Long, Y.-B. Wang, X.-F. Tan and W.-J. Zhang, Chem. Phys. Lett., 2013, 581, 26 CrossRef CAS; (d) B. Long, Z.-W. Long, Y.-B. Wang, X.-F. Tan, Y.-H. Han, C.-Y. Long, S.-J. Qin and W.-J. Zhang, ChemPhysChem, 2012, 13, 323 CrossRef CAS PubMed.
  393. Y. Kim, J. Am. Chem. Soc., 1996, 118, 1522 CrossRef CAS.
  394. R. L. Bell and T. N. Truong, J. Chem. Phys., 1994, 101, 10442 CrossRef CAS.
  395. R. C. De, M. Oliveira and G. F. Bauerfeldt, J. Phys. Chem. A, 2015, 119, 2802 CrossRef PubMed.
  396. I. M. Alecu and P. Marshall, J. Phys. Chem. A, 2014, 118, 11405 CrossRef CAS PubMed.
  397. M. Kumar, J. M. Anglada and J. S. Francisco, J. Phys. Chem. A, 2017, 121, 4318 CrossRef CAS PubMed.
  398. C. S. Tautermann, M. J. Loferer, A. F. Voegele and K. R. Liedl, J. Chem. Phys., 2004, 120, 11650 CrossRef CAS PubMed.
  399. A. F. Voegele, C. S. Tautermann, T. Loertingy and K. R. Liedl, Phys. Chem. Chem. Phys., 2003, 5, 487 RSC.
  400. T. Loerting, A. F. Voegele, C. S. Tautermann, K. R. Liedl, L. T. Molina and M. J. Molina, J. Geophys. Res.: Atmos., 2006, 111, D14307 CrossRef.
  401. U. V. Bhandarkar, M. T. Swihart, S. L. Girshick and U. R. Kortshagen, J. Phys. D: Appl. Phys., 2000, 33, 2731 CrossRef CAS.
  402. U. V. Bhandarkar, U. R. Kortshagen and S. L. Girshick, J. Phys. D: Appl. Phys., 2003, 36, 1399 CrossRef CAS.
  403. S. J. Warthesen and S. L. Girshick, Plasma Chem. Plasma Process., 2007, 27, 292 CrossRef CAS.
  404. P. Agarwal and S. L. Girshick, Plasma Sources Sci. Technol., 2012, 21, 055023 CrossRef.
  405. L. Ravi and S. L. Girshick, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 026408 CrossRef PubMed.
  406. P. Agarwal and S. L. Girshick, Plasma Chem. Plasma Process., 2014, 34, 489 CrossRef CAS.
  407. P. Seal and D. G. Truhlar, J. Am. Chem. Soc., 2014, 136, 2786 CrossRef CAS PubMed.
  408. P. Seal, J. Zheng and D. G. Truhlar, J. Phys. Chem. C, 2015, 119, 10085 CAS.
  409. J. L. Bao, P. Seal and D. G. Truhlar, Phys. Chem. Chem. Phys., 2015, 17, 15928 RSC.
  410. I. Oueslati, B. Kerkeni, A. Spielfiedel, W.-U. L. Tchang-Brillet and N. Feautrier, J. Phys. Chem. A, 2014, 118, 791 CrossRef CAS PubMed.
  411. T. V. Albu, B. J. Lynch, D. G. Truhlar, A. C. Goren, D. A. Hrovat, W. T. Borden and R. A. Moss, J. Phys. Chem. A, 2002, 106, 5323 CrossRef CAS.
  412. T. Takayanagi and S. Koido, Comput. Theor. Chem., 2017, 1115, 4 CrossRef CAS.
  413. J. Ho, J. Zheng, R. Meana-Pañeda, D. G. Truhlar, E. J. Ko, G. P. Savage, C. M. Williams, M. L. Coote and J. Tsanaktsidis, J. Org. Chem., 2013, 78, 6677 CrossRef CAS PubMed.
  414. A. Galano, L. Muñoz-Rugeles, J. R. Alvarez-Idaboy, J. L. Bao and D. G. Truhlar, J. Phys. Chem. A, 2016, 120, 4634 CrossRef CAS PubMed.
  415. J. Zheng, Y. Zhao and D. G. Truhlar, J. Chem. Theory Comput., 2009, 5, 808 CrossRef CAS PubMed.
  416. C. Alhambra, J. Gao, J. C. Corchado, J. Villa and D. G. Truhlar, J. Am. Chem. Soc., 1999, 121, 2253 CrossRef CAS.
  417. Q. Cui and M. Karplus, Adv. Protein Chem., 2003, 66, 315 CrossRef CAS PubMed.
  418. P. F. Faulder, G. Tresadern, K. K. Chohan, N. S. Scrutton, M. J. Sutcliffe, I. H. Hillier and N. A. Burton, J. Am. Chem. Soc., 2001, 123, 8604 CrossRef CAS PubMed.
  419. C. Alhambra, M. L. Sánchez, J. C. Corchado, J. Gao and D. G. Truhlar, Chem. Phys. Lett., 2001, 347, 512 CrossRef CAS ; corrected printing: 2002, 355, 388.
  420. A. Kohen and J. P. Klinman, J. Am. Chem. Soc., 2000, 122, 10738 CrossRef CAS.
  421. L. S. Devi-Kesavan and J. Gao, J. Am. Chem. Soc., 2003, 125, 1532 CrossRef CAS PubMed.
  422. M. Garcia-Viloca, D. G. Truhlar and J. Gao, Biochemistry, 2003, 42, 13558 CrossRef CAS PubMed.
  423. J. Pu, S. Ma, J. Gao and D. G. Truhlar, J. Phys. Chem. B, 2005, 109, 8551 CrossRef CAS PubMed.
  424. J. Pang, J. Pu, J. Gao, D. G. Truhlar and R. K. Allemann, J. Am. Chem. Soc., 2006, 128, 8015 CrossRef CAS PubMed.
  425. M. Garcia-Viloca, C. Alhambra, D. G. Truhlar and J. Gao, J. Am. Chem. Soc., 2002, 124, 7268 CrossRef CAS PubMed.
  426. M. Garcia-Viloca, C. Alhambra, D. G. Truhlar and J. Gao, J. Comput. Chem., 2003, 24, 177 CrossRef CAS PubMed.
  427. T. D. Poulsen, M. Garcia-Viloca, J. Gao and D. G. Truhlar, J. Phys. Chem. B, 2003, 107, 9567 CrossRef CAS.
  428. M. Roca, J. Andres, V. Moliner, I. Tunon and J. Bertran, J. Am. Chem. Soc., 2005, 127, 10648 CrossRef CAS PubMed.
  429. I. Feierberg, V. Luzhkov and J. Åqvist, J. Biol. Chem., 2000, 275, 22657 CrossRef CAS PubMed.
  430. A. Dybala-Defratyka, P. Paneth, R. Banerjee and D. G. Truhlar, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 10774 CrossRef CAS PubMed.
  431. R. Banerjee, D. G. Truhlar, A. Dybala-Defratyka and P. Paneth, Hydrogen Atom Transfers in B12 Enzymes, in Hydrogen-Transfer Reactions, ed. J. T. Hynes, J. P. Klinman, H.-H. Limbach and R. L. Schowen, Wiley-VCH, Weinheim, Germany, 2007, vol. 4, p. 1473 Search PubMed.
  432. A. Krzemińska, V. Moliner and K. Świderek, J. Am. Chem. Soc., 2016, 138, 16283 CrossRef PubMed.
  433. J. Javier Ruiz-Pernía, M. Garcia-Viloca, S. Bhattacharyya, J. Gao, D. G. Truhlar and I. Tuñón, J. Am. Chem. Soc., 2009, 131, 2687 CrossRef PubMed.
  434. I. Lans, J. Ramón Peregrina, M. Medina, M. Garcia-Viloca, À. González-Lafont and J. M. Lluch, J. Phys. Chem. B, 2010, 114, 3368 CrossRef CAS PubMed.
  435. J. Pang, S. Hay, N. S. Scrutton and M. J. Sutcliffe, J. Am. Chem. Soc., 2008, 130, 7092 CrossRef CAS PubMed.
  436. J. Gao, S. Ma, D. T. Major, K. Nam, J. Pu and D. G. Truhlar, Chem. Rev., 2006, 106, 3188 CrossRef CAS PubMed.
  437. D. Sicinska, D. G. Truhlar and P. Paneth, J. Am. Chem. Soc., 2005, 127, 5414 CrossRef CAS PubMed.
  438. M. J. Sutcliffe and N. S. Scrutton, Phys. Chem. Chem. Phys., 2006, 8, 4510 RSC.
  439. I. Tejero, M. Garcia-Viloca, A. González-Lafont, J. M. Lluch and D. M. York, J. Phys. Chem. B, 2006, 110, 24708 CrossRef CAS PubMed.
  440. S. Ferrer, I. Tuñón, S. Martí, V. Moliner, M. Garcia-Viloca, À. González-Lafont and J. M. Lluch, J. Am. Chem. Soc., 2006, 128, 16851 CrossRef CAS PubMed.
  441. N. Kanaan, S. Ferrer, S. Martí, M. Garcia-Viloca, A. Kohen and V. Moliner, J. Am. Chem. Soc., 2011, 133, 6692 CrossRef CAS PubMed.
  442. M. P. T. Duong and Y. Kim, J. Phys. Chem. A, 2010, 114, 3403 CrossRef CAS PubMed.
  443. B. K. Mai and Y. Kim, Chem. – Eur. J., 2014, 20, 6532 CrossRef CAS PubMed.
  444. J. Ángel Martínez-González, M. González, L. Masgrau and R. Martínez, ACS Catal., 2015, 5, 246 CrossRef.
  445. I. Tejero, N. González-García, A. González-Lafont and J. M. Lluch, J. Am. Chem. Soc., 2007, 129, 5846 CrossRef CAS PubMed.
  446. A. Thomas, D. Jourand, C. Bret, P. Amara and M. J. Field, J. Am. Chem. Soc., 1999, 121, 9693 CrossRef CAS.
  447. Q. Cui and M. Karplus, J. Am. Chem. Soc., 2002, 124, 3093 CrossRef CAS PubMed.
  448. J. W. Storer and K. N. Houk, J. Am. Chem. Soc., 1993, 115, 10426 CrossRef CAS.

This journal is © The Royal Society of Chemistry 2017