Quantum Tunneling Rates of Gas-Phase Reactions from On-the-Fly Instanton Calculations

: The instanton method obtains approximate tunneling rates from the minimum-action path (known as the instanton) linking reactants to the products at a given temperature. An e ﬃ cient way to ﬁ nd the instanton is to search for saddle-points on the ring-polymer potential surface, which is obtained by expressing the quantum Boltzmann operator as a discrete path-integral. Here we report a practical implementation of this ring-polymer form of instanton theory into the Molpro electronic-structure package, which allows the rates to be computed on-the-ﬂ y, without the need for a ﬁ tted analytic potential-energy surface. As a test case, we compute tunneling rates for the benchmark H + CH 4 reaction, showing how the e ﬃ ciency of the instanton method allows the user systematically to converge the tunneling rate with respect to the level of electronic-structure theory.

A ts u fficiently low temperatures, reactions with barriers are dominated by tunneling. 1 For hydrogen-and protontransfer reactions, the crossover temperatures (below which tunneling dominates) can by quite high, e.g., about 327 K in the gas-phase H + CH 4 reaction. Instanton rate theory has become a popular approach for calculating such tunneling effects 2−7 and can be derived from a semiclassical approximation to the exact rate constant. 8 The instanton describes the optimal tunneling pathway through the reaction barrier, and corresponds to an unstable periodic orbit on the inverted potential-energy surface (PES). 2 An efficient way to find the instanton is to search for saddle points on the ring-polymer potential surface, 9−11 which is obtained by expressing the Boltzmann operator as an integral over discretized Feynman paths. 12 This ring-polymer instanton method (also referred to as harmonic quantum transition-state theory 9 ) has been used to compute tunneling rates for a wide range of reactions. 9,13−19 A major advantage of the semiclassical instanton method is that it requires relatively few evaluations of the PES. Thus, although the dynamics and statistics are described less accurately than in exact quantum dynamics, 20 or in fully quantum statistical methods such as ring-polymer molecular dynamics (RPMD) 21−23 and the quantum instanton method, 24,25 it can be used with direct on-the-fly computation of the PES, using a high level of electronic-structure theory. This makes good practical sense, since often the errors in the PES dominate a rate calculation, making it more important to converge the electronic-structure calculation than to finesse the dynamics. Here we report an implementation of the ringpolymer instanton method 10 in the Molpro electronic-structure package, 26,27 which can be used to compute tunneling rates in gas-phase reactions using any electronic-structure ansatz available in the package. The restriction to gas-phase kinetics is made because instanton theory is not usually applicable to liquids, where the tunneling does not pass through a single barrier; such systems are better treated using RPMD. 10,21−23 Before applying instanton theory, it is usual to first locate the transition state and hence find the height and curvature of the barrier. An approximation to the crossover temperature is given by T c = ℏω b /2πk B where iω b is the imaginary vibrational wavenumber at the transition state. For sufficiently simple reaction pathways, it is only below this temperature that instantons exist and the rate formulas presented below are valid. However, note that for more complicated molecular processes, a significant amount of corner-cutting 3 can occur and the instanton pathway does not pass near a transition state 28 and can thus also exist at higher temperatures. 29 It is well-known that unlike RPMD or the quantum instanton method, the semiclassical instanton approach overestimates the rate just below the crossover temperature. A number of methods for removing this problem have been suggested, 30−34 but are not considered here.
Unlike earlier instanton approaches, 7 it is not necessary for the ring-polymer instanton approach to use a reduceddimensionality model for the reaction. We consider molecular systems with f degrees-of-freedom, equal to 3 times the number of atoms, in Cartesian coordinates. The reciprocal temperature is β =1/k B T, and we use the most abundant isotopic masses m j throughout. The discretized pathways can be described by ring polymers, and the instanton is located by optimizing a saddle point on the ring-polymer potential, where x 0,j ≡ x N,j , β N = β/N, and N is the number of system replicas known as ring-polymer beads. V(x 1 ,...,x f ) is the Born− Oppenheimer PES, which can be computed by any available electronic-structure method. The instanton rate is obtained when the results have converged with respect to N. Note that a number of ways of further increasing the efficiency have been suggested by discretizing the instanton more flexibly 16,35−37 but are not implemented here. Because the instanton folds back on itself, the computational effort can be halved by taking into account that half of the beads lie directly on top of the other half such that only N/2 electronic-structure calculations are necessary. 17 The algorithm implemented in Molpro works with a modified set of equations to optimize the half instanton, from which the full solution is easily obtained. The instanton geometry x̃is optimized from an initial guessed configuration, using a quasi-Newton optimization algorithm 38 with a Powell update for the ring-polymer Hessian. 39 Before the first step, the gradients and f × f Hessians are computed for each bead, but only the gradients need be recomputed after each optimization iteration. The Hessian of the half-ring-polymer is a × Nf Nf 1 2 1 2 banded matrix, which can be stored and diagonalized efficiently. The implementation uses the existing Molpro framework for computing first and second energy derivatives, which makes use of analytical derivatives wherever these are available, but otherwise performs finite numerical differentiation. More information about how to set up the instanton calculations can be found in the relevant section of the package manual, along with some examples. 27 The number of optimization steps can be reduced dramatically by using initial configurations close to the instanton geometry. One starts with an initial configuration in which the ring-polymer contains only a few beads, and is stretched over the transition state, 17 where q is the eigenvector corresponding to the imaginary mode at the transition state x ‡ and the variable Δ is chosen by trial and error such that the initial guess is good enough for the optimization to converge. The optimized configuration thus yields a few-bead approximation to the instanton, which is then interpolated onto a denser grid of beads, to give an initial guess for another optimization of the instanton geometry, which typically requires only a few iterations to optimize. Other techniques have also proven beneficial such as using optimized configurations from higher temperature instantons or with lower levels of electronic-structure theory as initial guesses. In any case, the results found are not dependent on the initial guess as long as it is in the region of convergence around the stationary point. For gas-phase reactions, this is usually easy to achieve. The total computational cost to locate the instanton and compute the rate at a given temperature is about N/2 times (in this case, 64 times) greater than that of a standard transition-state search and frequency calculation. As the program is written in a parallel manner, this cost can easily be shared over a number of processors, typically one per bead.
Our implementation compares the instanton rate with a variation of an Eyring TST calculation, for which the quantum partition functions are substituted by those that would be obtained by an N-bead ring polymer. This tends to the usual Eyring rate in the limit N → ∞, but we use the ring-polymer version to compare with the instanton rate, as a certain amount of cancellation of errors will occur which improves the convergence with respect to N. This version of Eyring's TST rate is defined as where the reactant partition function, Q r ,isdefined in terms of translational, rotational, and vibrational contributions from the separated reactant molecules in the usual way, and partition functions with the ‡ symbol refer to a ring polymer collapsed at the transition state configuration, The translational partition function describing a ring-polymer of total mass NM, Note that we have chosen to normalize these ring-polymer partition functions differently from the usual expressions. The rotational contribution to the partition function, Q rot ,isdefined using the classical expression as if all the beads in the ring polymer made up a Nf/3-atom "supermolecule" at reciprocal temperature β N . The moments-of-inertia tensor is given by where r⃗ i,a is the displacement from the center of mass of the ring-polymer of the ath atom (with mass m a ) of the ith replica, and  the 3 × 3 identity matrix. The ring-polymer rotational partition function is computed using the appropriate formula [Note that the symmetry number does not appear here but instead manifests itself in a number of identical transition states.]: where I B is the value of the nonzero eigenvalues of I (for linear configurations). The vibrational contribution of the ring polymer to the partition function is given by 40 where βω βω ℏ̃=ℏ sinh N k N k 1 2 1 2 and ω k are the vibrational wavenumbers at the transition state. As well as the imaginary mode, the f 0 modes corresponding to translational and rotational degrees of freedom are excluded from the product. There can be 5 or 6 of these depending on whether the transition state is linear or nonlinear.

Letter
The formula for the instanton rate can be written in a similar way to eq 3, 10 and instanton theory can thus be thought of as an extension of Eyring transition-state theory (TST) to the deep-tunneling regime: where the instanton action is Here, the translational and rotational partition functions are defined as in eq 4 and eq 6 using the ring-polymer instanton configuration, but there are two important differences in the vibrational contributions around the instanton. The full massweighted ring-polymer Hessian is recomputed after the optimization has finished and diagonalized to give the instanton frequencies η k . Whereas a ring polymer collapsed at the TS has only f 0 zero modes corresponding to translations and rotations, an instanton should also have one more zero mode corresponding to the permutation of the beads. 10 The integral around this cyclic permutational mode leads to a factor including the term The vibrational partition function is then defined in the following way: Note that this includes the magnitude of the imaginary frequency η 1 and that the eigenvector corresponding to this mode gives an approximation for the optimal reaction coordinate for an RPMD calculation. 10 The tunneling-factor is defined by  (12) which gives the factor by which the rate is increased upon including tunneling through the barrier. It is this result that is directly computed by the algorithm in the Molpro package. 26,27 κ tun is independent of the reactant partition function and is also much less sensitive to errors in the barrier height than the rate constant itself.
The gas-phase reaction H + CH 4 → H 2 +C H 3 is a wellstudied system, often used to compare and test new methodological developments in rate theory. 41−43 In this Letter, we report thermal rate coefficients for this reaction computed using the ring-polymer instanton approach based on energies obtained from on-the-fly coupled-cluster calculations.
Following Wu et al., 43 we employed the single and double excitation coupled cluster method (CCSD) with a spinrestricted Hartree−Fock reference wave function and spinsymmetry constraints on cluster amplitudes (RCCSD), 44 and, in some calculations, the standard perturbative correction for the effect of connected triple excitations (RCCSD(T)). 45 The cc-pVTZ basis set 46 was used, together with explicit twoparticle correlations through the F12a ansatz. 47 The sensitivity of the results to basis-set incompleteness was explored by carrying out calculations also using the cc-pVDZ basis. In all cases, N = 128 beads were used, which converged the results with respect to N to at least 2%.
The most important factor in the Eyring TST rate formula is the barrier height, whereas the tunneling effect depends on the shape of the barrier rather than its height. Table 1 shows predictions for the height and curvature of the barrier from the various levels of theory. These data suggest that the barrier height is being slightly overestimated by the more approximate methods, which would therefore predict the rate to be more than a factor of 10 too small at 200 K. The curvature is also seen to vary by about 5% and we will study how this affects the tunneling factor.
In Table 2, tunneling factors are presented at four different temperatures. The crossover temperature, T c , is predicted to be 327 K by RCCSD(T)-F12/cc-pVTZ. Instanton theory is only applicable below this temperature, and as the temperature drops, the tunneling factor increases. The tunneling factors evaluated at the different levels of theory vary by a considerable amount at low temperature, differing by more than a factor of 10 in some cases. Thus, even if we were somehow to correct for the error in the barrier height, this error would still remain. Note that, as expected, the variation correlates well with the curvature, in that high curvature leads to more tunneling. In fact, these large variations can be assigned to the negative exponent of the tunneling factor (rather than the prefactor). The values of S/ℏ−βV ‡ for the four methods at 200 K are 4.90, 5.52, 3.86, and 4.12 respectively.
As a final comparison, we study how the instanton pathways vary at the different levels of theory. Figure 1 shows the potential along the pathway as a function of cumulative massweighted path-length, This is the effective barrier through which the system must tunnel, and according to the Hamilton-Jacobi principle,   The Journal of Physical Chemistry Letters Letter completely defines the action, S, 35−37 and can thus be used to explain the variation in tunneling factors. For example, Figure 1 shows that RCCSD(T)-F12 yields a shorter, narrower barrier than RCCSD-F12 and hence smaller tunneling factors. Note that in this reaction, the instanton passes close to the saddlepoint, and there is relatively little corner-cutting; other reactions (e.g., those involving heavy-light-heavy combinations of atoms) could show much stronger dependencies of the instanton path on the level of electronic-structure theory. Table 3 compares the computed instanton rates [In order to take into account the degeneracy of the reaction due to the indistinguishability of the H atoms, the instanton rate formula, eq 8, has been multiplied by 4.] with CVT/μ-OMT, RPMD and multiconfigurational time-dependent Hartree (MCTDH) results, which were previously computed using analytic PESs, 43,48 which had both been fitted to RCCSD(T)/cc-pVTZ data points. [This is similar but not equivalent to our method which uses F12.] To help compare results, note that WWM 43 has a barrier height of 62.47 kJ/mol and frequency 1414 cm −1 , whereas the CBE PES 48 has a barrier height of 62.80 kJ/mol and frequency 1480 cm −1 . Nonetheless, the MCTDH results at 250 K on the CBE PES are only about 20% lower than those computed on the more accurate WWM surface. As we also expect the variation between the WWM PES and the on-the-fly calculations to be small, the variations between the rates presented in the table can mostly be assigned to the differences in dynamical methods.
Taking into account that it is well-understood that the instanton rates will be overestimated at 300 K, close to the crossover temperature, we see that there is only about a 30% error when compared to the benchmark MCTDH/WWM result at 250 K. This level of agreement has also been seen previously between MCTDH and instanton calculations carried out on the less accurate PJEG surface. 9 There is little difference between the instanton and CVT/μ-OMT rates at T = 250 K, as this particular reaction appears to follow a fairly simple reaction pathway without significant corner cutting, although the CVT/ μ-OMT results will probably deviate at lower temperatures when the tunneling factor becomes more sensitive to the pathway. 9 Bearing in mind that RPMD usually underestimates deep-tunneling rates by up to a factor of 2, 10 comparison between the instanton and RPMD results shows that there is little recrossing occurring. This therefore suggests that the instanton approach should be valid for this system and that the new result at 150 K should be an accurate estimate of the true rate.
It is worth pointing out that the accuracy of the instanton rates of Table 3, obtained using the RCCSD(T)-F12/cc-pVTZ method is better than that of the RPMD and MCTDH results of ref 53 and ref 9 carried out on the more approximate PJEG surface, which introduced about a factor of 4 error. This of course to be expected, since the tunneling depends exponentially upon the action of the instanton path, and thus on the shape of the PES in the vicinity of the barrier, as well as the barrier height itself. 54 For many reactions, for which accurate analytic potentials are not available, we expect that an on-the-fly implementation of instanton theory, such as the one implemented in Molpro that we have presented here, should give the most reliable estimates of the tunneling rate. (2) Miller, W. H. Semiclassical limit of quantum mechanical transition state theory for nonseparable systems. J. Chem. Phys. 1975Phys. , 62, 1899Phys. −1906 (

The Journal of Physical Chemistry Letters
Letter