Ab initio instanton rate theory made efficient using Gaussian process regression

Abstract

Ab initio instanton rate theory is a computational method for rigorously including tunnelling effects into the calculations of chemical reaction rates based on a potential-energy surface computed on the fly from electronic-structure theory. This approach is necessary to extend conventional transition-state theory into the deep-tunnelling regime, but it is also more computationally expensive as it requires many more ab initio calculations. We propose an approach which uses Gaussian process regression to fit the potential-energy surface locally around the dominant tunnelling pathway. The method can be converged to give the same result as from an on-the-fly ab initio instanton calculation but it requires far fewer electronic-structure calculations. This makes it a practical approach for obtaining accurate rate constants based on high-level electronic-structure methods. We show fast convergence to reproduce benchmark H + CH₄ results and evaluate new low-temperature rates of H + C₂H₆ in full dimensionality at a UCCSD(T)-F12b/cc-pVTZ-F12 level.

---

1 Introduction

Transition-state theory (TST) has surely become the most popular method for evaluating reaction rates in gas-phase chemistry.¹ It has achieved this status due to its simplicity and the fact that it can be evaluated with efficient computational algorithms. Two geometry optimisations are needed, for the reactant and transition states and two Hessian calculations, one at each stationary point. As only a small number of electronic-structure calculations are needed to evaluate the TST rate, expensive high-level ab initio methods can be used. This is necessary to achieve a good prediction, as small errors in the PES lead to exponential errors in the rate. TST however is based on classical dynamics and neglects important quantum effects such as tunnelling,² which can dominate the mechanism of certain chemical reactions of interest.^3–6

Ring-polymer instanton theory has proved itself to be a useful and accurate method for computing the rate of a chemical reaction dominated by tunnelling.⁷ The method is based on a first-principles derivation from the path-integral representation of the quantum rate^8–10 and can be thought of as a quantum-mechanical generalisation of TST. A ring-polymer discretisation of the path integral allows a simple optimisation algorithm to be used for locating the dominant tunnelling pathway, known as the “instanton”.^10–13 As with TST, it is possible to combine the instanton method with ab initio electronic-structure calculations to evaluate the potential-energy surface (PES) on the fly.^14–18 When compared with benchmark quantum dynamics approaches applied to polyatomic reactions, the instanton method typically gives low-temperature rates within about 20–30% of an exact calculation on the same PES.^15,19 This is, in many cases, less than the error in the rate which can be expected to result from the best achievable convergence of the electronic Schrödinger equation, implying that the accuracy of instanton theory itself is not the major issue.

The ab initio instanton method is very efficient when compared with other quantum dynamics approaches, including path-integral molecular dynamics or wave-function propagation. However, it remains considerably more computationally expensive than a TST calculation. The major reason for this expense is that energies, gradients and Hessians of the PES are required, not just at the transition state, but for each ring-polymer bead along the instanton, of which about 100 may be required. For high-accuracy electronic-structure methods, such as those provided by coupled-cluster theory, gradients and Hessians are typically evaluated using finite-differences, and can thus consume a lot of computational power. If the ring-polymer instanton method is to become widely applied in place of TST, the number of ab initio points will need to be reduced to bring the computational expense down, closer to that of a TST calculation.

It is important that high-quality electronic-structure calculations are employed as results can be strongly-dependent on the PES and give significant errors when using cheaper and less-accurate surfaces.^19,20 One suggestion for decreasing the computational effort required is to run the instanton calculation using a low-level surface and partially correct the result using a few high-level single-point calculations along the optimised pathway.^14,21,22 This approach (termed the ‘dual-level instanton approach’) certainly improves results, but it cannot always been relied upon as, in certain cases, the location of the instanton pathway may vary considerably depending on the quality of the PES. One can also use Taylor series expansions around the stationary points to obtain an approximate instanton solution analytically.^23–27 These approaches also have the potential to break down when the instanton pathway exhibits strong corner-cutting behaviour and deviates significantly from the transition state.

The procedure which has generally been followed for ring-polymer molecular dynamics rate theory^28,29 or wave-function propagation methods^30,31 has been to use an analytical function for the PES which is fitted to approximately reproduce ab initio points on the surface. In particular much attention has been given to water potentials,^32–34 on which instanton calculations have also been carried out for comparison with high-resolution spectroscopy.^35,36 Despite improvements and automation of this procedure, it remains a difficult task to fit a global potential, and is often based on tens of thousands of ab initio points,³⁷ computations which we wish to avoid.

The reason why these fitting procedures are typically difficult to carry out in practice is because a PES is a complex high-dimensional function. For many applications, including molecular dynamics or wave-function propagation, it is important to have a globally-accurate PES. In particular, if non-physical minima exist in the PES, the dynamics could be attracted there and give nonsensical results. Instanton theory has a particular advantage in that it only requires knowledge of a small region of the PES, located along a line representing the dominant tunnelling pathway. This implies that it might be possible to fit a locally-accurate surface around this small region in an efficient manner, as represented in Fig. 1. In this way we ensure that no extrapolation is used, but only interpolation, which is expected to be well behaved.

Fig. 1 The only areas of the PES which need to be accurately known are those around the instanton pathway or the reactant minimum (in order to obtain their partition functions). In this image, they are represented by the coloured areas, whereas those that are not built into the GPR are unshaded. The blue points represent the beads along the instanton path, while the black points represent the reactant and the transition state. Note that the tunnelling pathway cuts the corner to explore a space far from the transition state.

In this paper, we describe how we use Gaussian process regression (GPR)³⁸ to fit a local representation of the PES and thereby obtain the instanton rate using only a small number of ab initio calculations. By converging the rate with respect to the number of electronic-structure calculations, it is possible to obtain the same results as ab initio instanton theory, for a fraction of the cost. In this way, our GPR approach is almost as efficient as a TST calculation, but has the accuracy of a fully-converged ab initio instanton calculation. We are then able to take advantage of recent developments in high-accuracy electronic-structure methods,³⁹ which might otherwise be too expensive for an on-the-fly calculation. A similar combination of GPR and path-optimisation has been used successfully by the group of Jónsson.^40,41 A number of new developments are necessary for our implementation, as instanton theory also requires accurate knowledge of Hessians along the path, and because we apply the approach to gas-phase reactions, we must account for rotational invariance.

In the following, we describe the background theory as well as the particulars of our implementation of the approach. The results are then presented for two applications and the convergence properties discussed.

2 Theory

The results in this paper are computed by combining together a number of different approaches. Ring-polymer instanton theory is used to evaluate the rate based on a GPR fit to the PES, which has a training set composed of coupled-cluster electronic-structure calculations. It will be necessary to transform some data between different coordinate systems to use an appropriate set for each part of the calculation. The instanton equations are defined with Cartesians, as are the inputs and outputs of the electronic-structure calculations, but the GPR is best built using internal coordinates to ensure that it is rotationally invariant. In this way we formally make no further approximations to the instanton theory and also avoid having to construct a kinetic-energy operator in curvilinear coordinates.

2.1 Ring-polymer instanton theory

In the ring-polymer version of instanton theory,¹⁰ the dominant tunnelling pathway is represented by a path discretised into N segments. The points where the segments begin and end are given by Cartesian coordinates, x_i, called “beads”. Because the instanton pathway folds back on itself, only one half of the path need be specified.^12,13 A path defined by a set of N/2 beads, {x₁,…,x_N/2}, has the associated half-ring-polymer potentialwhere x_i,j is the Cartesian coordinate of the ith bead in the jth nuclear degree of freedom with associated mass m_j. The number of degrees of freedom is 3n, where n is the number of atoms. The spring constants are defined by the temperature, T, such that β_N = β/N and β = (k_BT)⁻¹.

The instanton configuration is defined as the saddle point of eqn (1) and in practical applications it can be located using quasi-Newton geometry optimisers.^11–13 These require gradients of the target function at each iteration but use update formulae to avoid recomputing the Hessians.⁴² The gradient of the ring-polymer potential depends on the gradients of the underlying PES at each bead geometry. In the on-the-fly implementation, these are obtained directly from an electronic-structure package, but here they are derivatives of the GPR fitted potential.

Once the instanton pathway is optimised, the theory accounts for fluctuations up to second order. Thus in order to evaluate the rate, we require Hessians of each bead. Again, these can be computed by an electronic-structure package or from the GPR. The calculation of a Hessian is usually carried out using second-order finite-differences and is therefore on the order of 3n-times more expensive than a gradient calculation.

Under the instanton approximation, the rate is given by

where the action is S = 2β_NℏU_N/2 and explicit expressions for the instanton vibrational, rotational and translational partition functions are given in ref. 10. The result should be converged with respect to the number of beads, N. Typically on the order of N = 100 beads are used to obtain a rate converged to two significant figures.

2.2 CCSD(T)-F12 theory

For electronic structures where the independent particle model is qualitatively correct, electronic energies computed at the basis set limit CCSD(T) level of theory are expected to be accurate to better than 1 kcal mol⁻¹ for reaction barriers, 0.1 pm for structures and 5 cm⁻¹ for harmonic vibrational wavenumbers.⁴³ Until relatively recently, the cost associated with using the large basis sets traditionally required to access the basis set limit has prevented this high level of theory from being routinely used in quantum dynamics simulations, which typically require many thousands of energy evaluations. With the maturation of modern F12 explicitly correlated theory,⁴⁴ near basis set limit CCSD(T) energies can now be computed using small (triple-zeta) orbital basis, at a cost only 15% larger than a traditional CCSD(T)/TZ calculation, and quantum dynamics studies can be performed using near basis set limit CCSD(T) Born–Oppenheimer potential-energy surfaces on a routine basis.

In F12 theory the standard manifold of correlating orbitals |ab〉 that parameterise two-body correlation functions in pair theories is supplemented with one geminal basis function per occupied orbital pair ij, chosen to directly model the Coulomb hole in the first-order pair correlation function

The correlation factor f(r₁₂) is chosen to be a linear combination of Gaussians⁴⁵ fit to an exponential function⁴⁶ with a length-scale of 1a₀, appropriate for valence electrons, and the many electron integrals that arise due to the explicit dependence on the interelectronic distance, r₁₂, and the presence of the strong orthogonality projector, Q̂, are decomposed into one- and two-electron components by inserting in approximate resolutions of the identity.⁴⁷ The coefficients t_ij^ab are optimised in the presence of fixed geminal contributions,⁴⁸ to reduce the geminal basis-set superposition error,⁴⁹ with coefficients chosen to satisfy the first-order singlet and triplet cusp conditions.⁵⁰ Small but numerically expensive geminal contributions to the energy Lagrangian function are neglected if they rank higher than third order in perturbation theory,⁵¹ resulting in the CCSD(T)(F12*) approximation.³⁹ In this work we use the Molpro electronic structure package⁵² and are restricted to using the slightly less accurate CCSD(T)-F12b approximation⁵³ where the geminal contributions from third order ring diagrams are also neglected. Nevertheless, the CCSD(T)-F12b energies computed in a TZ basis set are within 0.2 kJ mol⁻¹ per valence electron of the CCSD(T) basis set limit and retain the intrinsic accuracy of the wavefunction ansatz.⁵⁴

2.3 Gaussian process regression (GPR)

Gaussian process regression is a machine learning algorithm which can be used to efficiently generate complex hypersurfaces with limited data.³⁸ Recent work has applied this technique for the construction of potential-energy surfaces^55–57 and determining minimum energy paths^40,41 at a much lower computational cost. In this paper, a local representation of the PES is constructed around the instanton pathway and used to evaluate reaction rate constants.

Before carrying out the construction of a local PES with GPR, we first note that we have to utilise an internal coordinate system that accounts for rotational invariance. We define this internal coordinate system as q = q(x), where x is a set of Cartesian coordinates. This transformation to a rotationally invariant coordinate system is defined in Section 2.4.

In the simplest case, the training set is composed of known values of the potential, V(q_j), at the M reference points {q₁,…,q_M}. This defines the column vector y with elements

y_j = ε + V(q_j), (4)

where ε is an energy shift chosen such that the average of these elements is approximately zero. Noting that the derivative of a Gaussian process is also a Gaussian process,^38,40,41 it is also possible to include gradients and Hessians into the training set as described in ref. 40.

The potential for an unknown point q* can be predicted from GPR as

where k(q_i, q_j) is a covariance function for the prior. We chose a squared-exponential covariance function with length-scale γ and a prefactor f:The elements, w_j, of the vector, w, are determined by solving the linear equations

(K + σ²I)w = y, (7)

where the covariance matrix is defined by K_ij = k(q_i, q_j). By differentiating eqn (5), one obtains expressions for the gradient and Hessian of the PES. Because the covariance function is smooth, the PES is guaranteed to be differentiable to all orders. This is in contrast to the otherwise similar method of Shepard interpolation which produces surfaces which are not smooth enough to carry out instanton optimizations.¹⁹

σ is a noise term, which is introduced to avoid overfitting, and should be chosen to be the expected self-consistent error in the reference data. Together, f, σ and γ are known as hyperparameters. Their values can be optimised by maximising the log marginal likelihood,

Alternatively, one can also optimise the hyperparameters through the minimisation of errors by cross-validation.³⁸

The method above allows us to construct a local PES from a training set of M points. In our implementation, the general idea is as follows. Firstly, we construct an approximate PES with GPR using a small number of points, and then optimise the ring polymer based on this PES. After this, we refine the PES by adding new ab initio evaluations of points along the previously predicted ring-polymer configuration. Using the refined PES, we obtain a new ring-polymer configuration and then compare it with the previous one to check if the pathway has converged to the true instanton pathway. If this is not satisfied, the PES is refined again through the addition of more ab initio evaluations; this is continued iteratively until the convergence is achieved. The abovementioned scheme is further elaborated in Section 3.

The general scheme described above is similar to that done by the group of Jónsson,^40,41 wherein they obtain the minimum energy path using a GPR-aided nudged elastic band (NEB) method. This appears to have been highly successful, effectively reducing the number of ab initio evaluations required by an order of magnitude in comparison to a conventional NEB calculation. In this paper, we intend to emulate this drastic reduction in the computational effort for locating the instanton pathway and evaluating rates. As mentioned before, there are some differences in our implementation, such as the need for rotational invariance and accurate knowledge of the Hessians. We have found that the accuracy of the Hessians returned by GPR is significantly improved by explicitly providing Hessian data into the training set.

2.4 Non-redundant internal coordinate system

We would like to build the GPR representation of the PES using an internal coordinate system which is rotationally and translationally invariant. This is necessary as the relative rotational orientation of individual beads along the instanton pathway is not known a priori. However, we will need to be able to convert the information obtained from the GPR-based PES back into a Cartesian coordinate system in order to evaluate the instanton rate. Also note that the data available from electronic-structure packages are in Cartesian coordinates, which will need to be converted into internal coordinates in order to build the GPR-based PES.

Much recent work into machine-learning algorithms for describing intermolecular forces has further required permutational invariance.^58–61 Such advanced approaches could also be applied to our problem. However, as we only need to fit the potential locally, it is an unnecessary complication and thus we choose to neglect permutational symmetry. For our studies here, this is no inconvenience as we only need to compute the instanton rate for one of the equivalent reaction pathways and multiply the rate by the degeneracy.

A simple translationally and rotationally-invariant coordinate system for representing molecular geometries is provided by the n × n distance matrix,⁶² defined as

where r⃑_i is a three-dimensional vector of the Cartesian coordinates of atom i, such that r⃑₁ = (x₁, x₂, x₃), r⃑₂ = (x₄, x₅, x₆), etc.

Although it is possible to convert data from a Cartesian coordinate system into this set,⁶³ the back transformation is not well defined, as the internal coordinates are redundant. In order to obtain a non-redundant set of internal coordinates, we follow the approach of Baker et al.⁶⁴ Firstly, we unravel the matrix D to give the coordinates as a vector of length n²,

d = [D₁₁D₁₂ … D₂₁D₂₂ … D_nn]^T. (10)

The B matrix is defined to describe how changes in the Cartesian coordinates affect these redundant coordinates as

The elements of this n² × 3n matrix are given explicitly by

where α runs over the indices of three-dimensional space.

A square matrix, G = BB^T, is formed and then diagonalised to obtain the eigenvalues and eigenvectors. The non-redundant eigenvectors are those corresponding to the nonzero eigenvalues (of which there will be 3n − 6 for a nonlinear isolated molecule), whereas the redundant eigenvectors have zero eigenvalues. The non-redundant eigenvectors are collected into the columns of a matrix, U. With this, we can now transform d into a non-redundant coordinate system, defined by

q = U^Td. (13)

It is this internal coordinate system which is used to build the GPR representation.

Note that the matrix U is built only once at a reference geometry and is used to define the transformation to q at all other geometries. The reference geometry used in our studies was the transition state, although this is not a requirement. The same U matrix is then used for new geometries to give a consistent definition of the internal coordinates q = q(x).

Therefore the required relationship between the internal coordinates and Cartesians is given by dq = B_qdx, where

The gradient and Hessian in the non-redundant internal coordinate system are defined as

and similarly, in Cartesian coordinates,

Given a geometry x to define the appropriate orientation, the gradients and Hessians obtained from the GPR in internal coordinates can be transformed back to Cartesian coordinates. Obtained using the chain rule, the transformations are defined by

g_x = B_q^T g_q (17)

where .

In order to transform the gradients and Hessians obtained from electronic-structure calculations into the q coordinate system, these equations need to be inverted. However, as B_q is not a square matrix, we need to define the generalised inverse as

(B_q^T)⁻¹ = (B_qB_q^T)⁻¹B_q. (19)

The required transformations are

g_q = (B_q^T)⁻¹g_x (20)

These equations define all of the necessary transformations needed for converting the ab initio data into reduced coordinates, and for converting it back to a Cartesian system at a given orientation.

3 Method

Our aim is to reproduce the same result as an ab initio instanton calculation performed on the fly. As with these calculations, we must therefore consider convergence with respect to N. For our new approach based on GPR, we must also simultaneously converge the result with respect to the number of points in the training set.

Here we outline our standard protocol for computing converged instanton rates using GPR. This is made up of two parts: first, in which the instanton pathway is located, and second, in which the fluctuation terms are converged to yield the final instanton rate. We have attempted to design this protocol to be stable and efficient. In our study, we have found that this protocol posed no significant problems for the systems tested here. In future studies, one could consider improvements which may increase the efficiency further. In a realistic working environment, a researcher has the freedom to add information to the GPR however they like until the result is converged.

Our protocol is designed for the case that single-point ab initio calculations are by far the most expensive part of the calculation. We also assume that the Hessian calculations are orders of magnitude slower than potential or gradient evaluations. This is commonly the case for many electronic-structure methods, especially if the Hessians are computed using finite differences. The efficiency of our protocol should thus be measured in terms of the number of ab initio calculations required, and in particular the number of Hessians. We present these numbers for specific examples in the next section.

The protocol described below is intended for a calculation of a single instanton rate at a given temperature, as is the approach used in the H + CH₄ benchmarks we present in Section 4.1. If, as is common, one needs the rate at multiple temperatures, it is recommended to start just below the crossover temperature, T_c. The optimised instanton can be used as the initial guess and GPR training set for a calculation at a lower temperature. We use this more efficient approach for our H + C₂H₆ calculations in Section 4.2.

3.1 Protocol

(1) Optimise the reactants and transition state (using a standard Quantum Chemistry package), and obtain gradients and Hessians for the optimised geometries. The optimised transition-state geometry in Cartesian coordinates is notated x^‡.

(2) By diagonalising the mass-weighted Hessian at the transition-state, calculate the cross-over temperature,

where ω_b is the magnitude of the imaginary frequency.

(3) An initial guess for the instanton configuration is obtained using¹³

where z is the normalised non-mass-weighted eigenvector corresponding to the imaginary mode at the transition state and Δ is a user-defined spread of points. Typically we choose Δ ∼ 0.1 Å and N = 16 for an initial guess.

However, if previous instanton optimisations at a higher temperature have been performed successfully, these configurations usually provide a better initial guess.

(4) Calculate ab initio potentials and gradients for the points obtained in step 23.

(5) Repeat until convergence:

(a) Optimise hyperparameters using methods defined previously under Section 2.3.

(b) Starting with a low number of beads N, locate the ring-polymer path by increasing N until the action S/ℏ is converged to 2 decimal places.

(c) Check if the mean bead displacement , where PC corresponds to the path convergence limit. Also check that the convergence of the action |S^new − S^old|/ℏ ≤ 10⁻².

• If this is satisfied, this means that the ring-polymer path has converged. Continue to step 6.

• Otherwise provide new inputs (ab initio energies and gradients) to the GPR training set along the current ring polymer (i.e. increase the number of training points M) and then go back to step 5a.

(6) Repeat until convergence is achieved:

(a) Provide a couple of new points along the converged instanton pathway to the GPR training set, this time also including the Hessians.

(b) Optimise the hyperparameters using the methods defined previously under Section 2.3.

(c) Locate the instanton pathway and calculate the rate, k, increasing N until this converges.

(d) Test if |k^new − k^old|/k^new ≤ RC, where RC corresponds to the rate convergence limit.

• If this is satisfied, the iterative algorithm is terminated, and the current value of k is taken as the converged instanton rate.

• Otherwise, return to step 6a.

In the following calculations, we built the GPR using energies in hartrees (E_h) and Cartesian coordinates in ångströms (Å). In these units, the typical values used for the length-scale were γ ∼ 0.3–0.4, and for the prefactor, f = 0.09. We specified the noise term differently for the potentials, gradients and Hessians, as σ_V ∼ 10⁻⁶, σ_G ∼ 10⁻⁴ and σ_H ∼ 10⁻³. The convergence limits of PC = 10⁻² Å and RC = 10⁻² were used.

We have outlined the simplest protocol which has the desired properties of converging the instanton rate without needing a large number of ab initio calculations. However, it is not necessarily the optimal choice for all problems. In particular, it should be noted that in this work, new information is provided to the GPR training set at the positions of beads chosen by hand. This was done in a systematic way, wherein during the path convergence step, the beads were chosen such that they are evenly distributed along the current ring polymer. Once the path is converged, beads where the Hessians are to be included were chosen in a similar manner, i.e. evenly distributed along the converged pathway. There may be better ways of providing new information to the GPR training data; for instance one can evaluate the expected fitting error along the current pathway and then provide points at the areas with high variance. By being more selective, one can potentially further reduce the number of ab initio calculations required.

4 Results

The method described above was applied to the following two systems:

H + CH₄ → H₂ + CH₃

H + C₂H₆ → H₂ + C₂H₅.

The first is a standard benchmark reaction for testing quantum rate theories and has been studied with various methods including MCTDH,^65,66 ring-polymer molecular dynamics,⁶⁷ the quantum instanton,⁶⁸ as well as ring-polymer instanton theory.^12,15,19

The second reaction is beyond the current limits of exact quantum mechanics unless reduced dimensionality models are used. Using the GPR formalism, we are able to present a converged ab initio instanton rate for the first time. We compare these results with those predicted by other semiclassical methods.

4.1 H + CH₄

An on-the-fly ab initio instanton calculation has been done by one of us for this polyatomic reactive system.¹⁵ Here, we use this reaction as a benchmark case for our GPR-aided instanton calculation and show that we are able to obtain the same result as an on-the-fly calculation with a significant reduction in the number of potentials, gradients, and most importantly, Hessians required.

The electronic-structure method used in ref. 15 was RCCSD(T)-F12a/cc-pVTZ and we use exactly the same method for the training set for the GPR. Note that in this paper, as well as in ref. 15, this method is also used to obtain the properties of the isolated reactants (including the H atom). To account for the indistinguishability of the H atoms, the instanton rate formula is multiplied by 4.

The results in ref. 15 were computed using N = 128, which we also use here. This required the calculation of 64 ab initio potentials and gradients per iteration of the instanton optimisation scheme. Because approximately 10 iteration steps are usually required for an instanton optimisation, about 640 gradients were computed in addition to the 64 Hessians once the instanton had been optimised. Here, we followed the protocol outlined in the previous section independently for three different temperatures. This allows us to accurately determine the computational effort required for a converged rate.

In Table 1, the rows correspond to the iterations of step 5 of the protocol, in which the pathway is optimised by adding more potentials and gradients to the GPR training set. The action is seen to converge to two decimal places after only a few iterations. Here, this was done with fewer than 50 potentials and gradients for all three temperatures. This means that a reduction in the number of gradient evaluations by an order of magnitude has been achieved.

Table 1 Convergence of the instanton path with the iteration of protocol step 5. The number of potentials (V), gradients (G), and Hessians (H) included in the GPR training set is explicitly noted. Here the single Hessian in the training set corresponds to that of the transition state

T (K)	Iteration	Training set	S/ℏ	Δx (10^−3 Å)
300	1	10V, 10G, 1H	25.167	—
	2	19V, 19G, 1H	25.243	41.7
	3	28V, 28G, 1H	25.249	1.78
250	1	10V, 10G, 1H	28.957	—
	2	19V, 19G, 1H	29.244	37.3
	3	28V, 28G, 1H	29.301	3.17
	4	37V, 37G, 1H	29.274	1.77
	5	46V, 46G, 1H	29.278	0.13
200	1	10V, 10G, 1H	32.586	—
	2	19V, 19G, 1H	33.871	70.1
	3	28V, 28G, 1H	33.945	2.54
	4	37V, 37G, 1H	33.904	1.90
	5	46V, 46G, 1H	33.907	0.70

---

This fast convergence is also observed in Fig. 2, where it is seen that, at the lowest temperature studied, the pathway already has the correct shape after the second iteration. In this figure, the potential along the pathway is plotted as a function of a cumulative mass-weighted path length,

Fig. 2 Convergence of the ring-polymer instanton at 200 K for H + CH₄. The initial GPR training set was defined by eqn (23). The ring-polymer beads are plotted as a function of their potential energy and the path length, l, as defined by eqn (24).

It should be noted that the plots are shifted such that they are centred around l = 0.

In Table 2, the GPR model is further refined, as described in step 6 of the protocol, by providing more observations (i.e. more ab initio potentials, gradients and Hessians) to the GPR training set. Our findings show that it is necessary to include a few Hessians directly in the GPR training set, and that the transition-state Hessian alone is not sufficient to describe the fluctuation terms of the instanton. Note that at low temperatures, the GPR requires a few more Hessians to converge the rate. This is due to the fact that the instanton stretches out more at lower temperatures, thus meaning that the GPR needs more information as the instanton covers a larger area of the PES.

Table 2 The rates obtained from the GPR-based instanton calculations are given as the information provided to the GPR training set is increased. The error is measured relative to the on-the-fly ab initio results of ref. 15. Note that for the rate calculation, one further Hessian is needed at the reactant geometry, but that this is not included in the GPR training set

T (K)	Training set	k (cm^3 s^−1)	Relative error
300	28V, 28G, 1H	5.49(−19)	220%
	31V, 31G, 4H	1.72(−19)	1.2%
	33V, 33G, 6H	1.69(−19)	<1%
	Ref. 15	1.70(−19)	—
250	46V, 46G, 1H	4.25(−20)	790%
	49V, 49G, 4H	4.74(−21)	−1.3%
	51V, 51G, 6H	4.80(−21)	<1%
	Ref. 15	4.80(−21)	—
200	46V, 46G, 1H	1.68(−20)	>1000%
	49V, 49G, 4H	0.98(−22)	−10%
	51V, 51G, 6H	1.07(−22)	−1.8%
	53V, 53G, 8H	1.08(−22)	<1%
	Ref. 15	1.09(−22)	—

---

The convergence is fast and it takes no more than 8 Hessians to converge the rates for all of the temperatures to less than 1% of that of the ab initio calculation. This is a remarkable improvement in terms of the computational effort required over the ab initio instanton calculations as the Hessian calculations account for a huge percentage of the computational effort required. Having reduced the number of Hessians required from 64 to 8, the reduction in computational power needed would allow us to investigate problems involving larger molecules and to also use higher-level electronic-structure methods.

4.2 H + C₂H₆

The H abstraction reaction from ethane follows the same mechanism as abstraction from methane. From a theoretical point of view, it is of interest as the number of degrees of freedom is significantly higher such that full-dimensional exact quantum methods are not applicable and approximations must be made. There are two types of approximations which can be used to make the simulation tractable. One makes use of semiclassical dynamics, and the other involves reducing the dimensionality of the system. The instanton method is an example of the former, as are other semiclassical extensions of the transition-state theory^25,69 and ring-polymer molecular dynamics.⁷⁰ Reduced-dimensionality models allow the quantum scattering theory to be applied⁷¹ and can also be combined with semiclassical approaches.^25,26,72 Experimental results are available at 300 K,^73,74 but unfortunately not at lower temperatures, where the tunnelling effect is more important. Here, we compare the results of our instanton rate calculations with other theoretical calculations, and discuss the relative efficiency of the various methods.

4.2.1 Ab initio calculations. Due to the efficiency of the GPR-aided instanton approach seen in our benchmark tests, we are able to use high-accuracy and computationally expensive electronic-structure methods. The method we choose is UCCSD(T)-F12b as discussed in Section 2.2. Table 3 shows the predictions for barrier heights, V^‡, and imaginary frequencies, ω_b, with increasingly large basis sets. Hessians with cc-pVQZ and cc-pV5Z basis sets were not evaluated due to the large amount of computational resources that would be required. However, we can see that cc-pVTZ-F12 reproduces almost the same barrier height as cc-pV5Z, in accordance with the study by Spackman et al.⁷⁵ which suggested that the cc-pVnZ-F12 basis sets have similar performance to the cc-pV(n + 2)Z basis sets (where n = D, T, Q, etc.) in terms of results when using CCSD(T)-F12. Hence in the following calculations, we will use the cc-pVTZ-F12 basis set.

Table 3 Barrier heights and imaginary frequencies for H + C₂H₆ using increasingly larger basis sets at the UCCSD(T)-F12b level

Method	V ^‡ (kJ mol^−1)	ω b (cm^−1)
UCCSD(T)-F12b/cc-pVDZ	54.57	1398
UCCSD(T)-F12b/cc-pVTZ	51.24	1461
UCCSD(T)-F12b/cc-pVTZ-F12	50.03	1469
UCCSD(T)-F12b/cc-pVQZ	50.46	—
UCCSD(T)-F12b/cc-pV5Z	50.07	—

---

With our chosen method, the crossover temperature is predicted to be 337 K. We ran three instanton calculations, first at 300 K, and then used this as a starting point for a calculation at 250 K, and in the same way for 200 K. This approach may slightly reduce the number of iterations needed for convergence. For instance, it can be seen in Fig. 3 that the optimisation of the 200 K instanton is obtained in only a few iterations and that the path is almost correct even after the first. The convergence criteria used for this system were similar to those used in the H + CH₄ system.

Fig. 3 Convergence of the ring-polymer instanton at 200 K for H + C₂H₆. The initial GPR training set was given by points along the 250 K instanton path. The path length, l, is defined by eqn (24).

The Cartesian representation of the optimised path for H abstraction from ethane is shown in Fig. 4. It is seen that the mechanism is similar to that of H + CH₄, shown in ref. 15, in that the abstracted hydrogen does most of the tunnelling, and is accompanied by a small movement of its neighbouring hydrogens. The atoms at the far end of the ethane molecule hardly participate in the instanton at all. Note, however, that they still make a contribution to the fluctuations, and thus cannot be neglected.⁷⁶

Fig. 4 Representation of the ring-polymer instanton for H + C₂H₆ at 200 K.

The results of our GPR-based instanton calculations are presented in Table 4. These rates account for the degeneracy of the reaction by multiplying the formula in eqn (2) by a factor of 6. The 300 K result was obtained with a training set including 33 potentials and gradients, and 6 Hessians. The calculations at the lower temperatures of 250 K and 200 K added an additional 6 Hessians to the training set (i.e. at 250 K, the training data includes 6 Hessians from 300 K and 6 Hessians from 250 K) in order to converge the rates. This represents a reduction in the computational effort by an order of magnitude, similar to what has been observed for H + CH₄.

Table 4 Calculated rates (in cm³ s⁻¹) for H + C₂H₆ obtained by the GPR-aided instanton method and other direct dynamics methods. The tunnelling factor, κ_tun, is defined as the ratio between the instanton rate and Eyring TST

T/K	GPR-aided instanton		SCTST^25	RD-QS^71
	κ tun	Rate
300	15	7.0(−17)	3.88(−17)	6.23(−17)
250	38	6.4(−18)	9.51(−18)	7.97(−18)
200	623	5.7(−19)	2.50(−19)	6.69(−19)

---

It is clear from our calculation that the tunnelling effect makes a large contribution to the rate, even at 300 K. This is confirmed by experimental results at this temperature, which in various setups have been measured to be 3.13 × 10⁻¹⁷ cm³ s⁻¹ (ref. 73) or 7.47 × 10⁻¹⁷ cm³ s⁻¹,⁷⁴ and which both lie in the same order of magnitude as our prediction. Note that we expect the instanton approach to slightly overpredict the rate (by up to a factor of 2) at 300 K as this lies close to the value of T_c.⁷⁷ Unfortunately, no experimental results are available for comparison at lower temperatures where the tunnelling effect is predicted to increase dramatically.

Table 4 also compares our predicted rate with those of the reduced-dimensionality quantum scattering (RD-QS) calculations by Horsten et al.⁷¹ and a full-dimensional semiclassical transition-state theory (SCTST) rate calculation by Greene et al.²⁵

The RD-QS calculations utilised a similar electronic structure method to our calculations, albeit with F12a rather than F12b, which gives a barrier height only 0.1 kJ mol⁻¹ lower. The SCTST calculations employed the CCSD(T)/cc-pVTZ method for the energies at the stationary points, which gives a barrier height 0.4 kJ mol⁻¹ lower. We expect these differences to lead to only a minor deviation.

The instanton results are in quite close agreement with RD-QS, where the rates differ by no more than 25%. This is what is typically expected when comparing results obtained with the instanton method and those obtained with exact quantum methods.¹⁹ This confirms that, at least for this system, the reduced-dimensionality approach is not causing an appreciable error in the tunnelling effect.

There is a slightly larger discrepancy between the instanton and SCTST results,²⁵ which increases at lower temperatures. The SCTST rate calculation involved a total of 118 ab initio Hessians at the MP2/cc-pVTZ level, with energies at stationary points evaluated with CCSD(T)/cc-pVTZ. The GPR-instanton method required only 6 Hessians to converge the rate at each temperature and thus a high-level of theory for the Hessian calculations can be used as well. There are two reasons for the discrepancy in the SCTST rates. One is that lower-level electronic-structure theory was used in ref. 25 for the Hessian calculations. The second is that at low temperatures, the instanton pathway stretches far from the transition state and the PES cannot therefore be well represented by a Taylor series around the transition state.

In this case, there are no dramatic differences between the theoretical predictions. It seems that the H + C₂H₆ reaction follows a simple pathway for which reduced-dimensionality models are applicable. However, we expect that for more complex reactions there will be a larger discrepancy and that, in many cases, the full-dimensional instanton theory will be the most accurate.

4.2.2 Results on a fitted PES. In order to get an idea of the accuracy of the instanton approach for this reaction, we compare instanton rates with those of other semiclassical approaches based on the fitted, global CVBMM potential-energy surface.⁶⁹

This PES was constructed by dividing the system into a reactive part which would be treated with semiempirical valence bond theory and a non-reactive part treated with molecular mechanics. It was parameterised against density functional theory, of which more details can be found in ref. 69. The barrier height obtained with the CVBMM PES is 47.90 kJ mol⁻¹ and has a predicted crossover temperature of 352 K.

Table 5 presents the rates of three methods, the instanton theory (this work), the quantum instanton theory (QI)⁷⁸ and the small curvature tunnelling correction to canonical variational TST (CVT/SCT).⁶⁹ The tunnelling factors are seen to be about a factor of 2 larger than those from the ab initio method, mainly due to the fact that the CVBMM barrier is too narrow and thus overpredicts the tunnelling factors

Table 5 Rate comparison between different methods using the CVBMM PES. All of the rates are in cm³ s⁻¹

T/K	Instanton		CVT/SCT^69	QI^78
	κ tun	Rate
300	24	1.25(−16)	1.44(−16)	1.15(−16)
250	80	1.46(−17)	—	1.40(−17)
200	1296	1.61(−18)	1.90(−18)	1.16(−18)

---

The CVT/SCT rate is in close agreement with that of the instanton theory, which at least in this cases implies that the dominant tunnelling pathway is well approximated by the minimum-energy pathway used by CVT/SCT. It is expected that, in general for more complex reactions, the instanton method, which defines the tunnelling pathway in a rigorous manner, will give a more accurate result.

Unlike the ring-polymer instanton approach, the QI method does not use a steepest-descent approximation and thus includes anharmonic vibrational effects in full dimensionality. In order to do this, it samples over a statistically large number of path-integral configurations and would therefore not be a practical computational method when combined with high-level ab initio potentials. Nonetheless, these anharmonic effects only change the rate by less than 50% at the lowest temperature studied. This is in agreement with the findings of ref. 78 which showed that, at low temperatures, a small increase in the rate resulted from making a harmonic approximation to the internal rotation. This confirms that instanton theory gives a reliable prediction of the order-of-magnitude of the rate. The real advantage of the instanton approach over this method is that it can be applied to new reactions without needing to build a global PES at all.

5 Conclusions

We have demonstrated how ab initio instanton theory can be made efficient using GPR to fit the PES locally around the dominant tunnelling path. This was demonstrated first using the H + CH₄ reaction as a benchmark, for which we have shown that the number of electronic-structure calculations can be reduced by an order of magnitude, while converging the rate to within 1% of the benchmark result. We then proceeded to evaluate the instanton rates for H + C₂H₆, based on UCCSD(T)-F12b/cc-pVTZ-F12 electronic-structure calculations. Most importantly, the number of Hessians needed for all of these calculations is about 6, which makes the method more efficient than full-dimensional SCTST calculations and almost as efficient as a classical TST calculation.

When studying a complex network of reactions, TST is commonly used to obtain a rate for the many possible reaction steps.⁷⁹ By evaluating the crossover temperature for each step, it can be easily determined whether tunnelling is likely to play a role, and instanton calculations can be run for these steps only. As there are typically many more steps for which tunnelling is not important than those for which it is, the number of ab initio calculations needed for the instanton calculations would be small in comparison to the overall total. In this way, tunnelling can be rigorously accounted for without significantly increasing the computational effort.

In this work, we suggested a simple protocol which, in our tests, showed no particular problems. We note, however, that it could still be improved in a number of ways which would further increase the efficiency. For instance, when using estimates of the GPR fitting error, we could select new points to be added to the training set in a more systematic way. These could also be used to estimate the fitting error in the rate constant in a similar way to what has been done for TST calculations.⁸⁰

Other techniques might allow us to reduce the number of high-level calculations by including low-level ab initio information into the GPR training set. One possibility would be to use this low-level information only for the initial iterations to locate the region of space where the instanton is likely to exist on the high-level surface. The final iteration could be done using only the high-level information to ensure convergence to the correct result. However, one could also consider combining the high- and low-level information in the training set, as in the dual-level approach.²² By using a larger value of the noise term for the low-level points, the GPR would then fit itself accurately to the high-level points, and the low-level information could be used as a rough guide for the shape. Typically the frequency calculations from the low-level calculations are a good approximation even if the absolute energies are not, and so most Hessians could be derived from the low-level calculations. One could imagine systematically converging to the correct result by adding more high-level ab initio points such that the accuracy would not be compromised.

We have shown in this paper that we can converge the rate with respect to the number of ring-polymer beads, N, as well as with respect to the number of points included in the GPR training set. However, the accuracy of our method is still limited by the computational expense of electronic-structure methods, which are rarely possible to fully converge. Methods such as F12 have been very useful for increasing this efficiency⁴⁴ as well as linear-scaling methods⁸¹ and the use of graphical processing units.⁸² Our GPR-aided instanton theory can take advantage of such improvements to ab initio electronic-structure calculations in an efficient way.

We did not find particularly large differences in the rate predictions for the H + C₂H₆ reaction between the instanton approach and other theories. This is due to the rather simple mechanism exhibited by the H abstraction reaction, which follows a pathway close to the minimum-energy path, making the CVT/SCT and reduced-dimensionality models valid. The advantage of the instanton theory is that no a priori choice of reduced coordinates, or tunnelling coordinate is made. This makes the approach applicable also to more complex reactions as well as tunnelling splitting calculations.³⁵ In these cases it is expected that the instanton path will deviate more strongly from the minimum-energy path, and the full-dimensional instanton theory will be required to obtain an accurate prediction. The proof of principle outlined in this work for combining GPR with the instanton theory will then be exploited in future studies of new reactions.

Note added in proof

It has come to our attention that in a recent paper, Kästner and coworkers apply a similar scheme to compute ab initio instanton rates using a neural network fit of the PES.⁸³ Their findings are similar to our own, in that a huge computational saving results from the approach and they also stress the importance of including ab initio Hessians in the training data. Neural networks have been found to be one of the best methods for fitting global PESs. However, our simpler GPR scheme may have some advantages over the neural networks in this case where only a local description is required based on a small training set. In particular, GPR has a limited number of hyperparameters which can be easily optimised, whereas neural networks appear to have a strong dependence on randomly chosen initial weights and require an averaging procedure to obtain a prediction for the rate.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been financially supported by the Swiss National Science Foundation (project no. 175696).

Notes and references

1. D. G. Truhlar, B. C. Garrett and S. J. Klippenstein, J. Phys. Chem., 1996, 100, 12771–12800.
2. E. Wigner, Trans. Faraday Soc., 1938, 34, 29–41.
3. R. P. Bell, The Tunnel Effect in Chemistry, Chapman and Hall, London, 1980.
4. B. K. Carpenter, Science, 2011, 332, 1269.
5. D. Ley, D. Gerbig and P. R. Schreiner, Org. Biomol. Chem., 2012, 10, 3781–3790.
6. J. Meisner and J. Kästner, Angew. Chem., Int. Ed., 2016, 55, 5400–5413.
7. J. O. Richardson, J. Chem. Phys., 2018, 148, 200901.
8. W. H. Miller, J. Chem. Phys., 1975, 62, 1899–1906.
9. J. O. Richardson, J. Chem. Phys., 2016, 144, 114106.
10. J. O. Richardson, Int. Rev. Phys. Chem., 2018, 37, 171.
11. J. O. Richardson and S. C. Althorpe, J. Chem. Phys., 2009, 131, 214106.
12. S. Andersson, G. Nyman, A. Arnaldsson, U. Manthe and H. Jónsson, J. Phys. Chem. A, 2009, 113, 4468–4478.
13. J. B. Rommel, T. P. M. Goumans and J. Kästner, J. Chem. Theory Comput., 2011, 7, 690–698.
14. G. V. Mil’nikov, K. Yagi, T. Taketsugu, H. Nakamura and K. Hirao, J. Chem. Phys., 2004, 120, 5036–5045.
15. A. N. Beyer, J. O. Richardson, P. J. Knowles, J. Rommel and S. C. Althorpe, J. Phys. Chem. Lett., 2016, 7, 4374–4379.
16. V. Ásgeirsson, A. Arnaldsson and H. Jónsson, J. Chem. Phys., 2018, 148, 102334.
17. T. P. M. Goumans and J. Kästner, Angew. Chem., Int. Ed., 2010, 49, 7350–7352.
18. M. Kryvohuz, J. Chem. Phys., 2012, 137, 234304.
19. K. Karandashev, Z.-H. Xu, M. Meuwly, J. Vaníček and J. O. Richardson, Struct. Dyn., 2017, 4, 061501.
20. K. Yagi, G. V. Mil’nikov, T. Taketsugu, K. Hirao and H. Nakamura, Chem. Phys. Lett., 2004, 397, 435–440.
21. G. V. Mil’nikov, K. Yagi, T. Taketsugu, H. Nakamura and K. Hirao, J. Chem. Phys., 2003, 119, 10–13.
22. J. Meisner and J. Kästner, J. Chem. Theory Comput., 2018, 14, 1865–1872.
23. W. H. Miller, R. Hernandez, N. C. Handy, D. Jayatilaka and A. Willetts, Chem. Phys. Lett., 1990, 172, 62–68.
24. T. L. Nguyen, J. F. Stanton and J. R. Barker, Chem. Phys. Lett., 2010, 499, 9–15.
25. S. M. Greene, X. Shan and D. C. Clary, J. Chem. Phys., 2016, 144, 244116.
26. S. M. Greene, X. Shan and D. C. Clary, J. Chem. Phys., 2016, 144, 084113.
27. Z. Smedarchina, W. Siebrand and A. Fernández-Ramos, J. Chem. Phys., 2012, 137, 224105.
28. R. Collepardo-Guevara, Y. V. Suleimanov and D. E. Manolopoulos, J. Chem. Phys., 2009, 130, 174713.
29. Y. V. Suleimanov, F. J. Aoiz and H. Guo, J. Phys. Chem. A, 2016, 120, 8488–8502.
30. S. C. Althorpe and D. C. Clary, in Annu. Rev. Phys. Chem., ed. S. R. Leone, Annual Reviews, Palo Alto, Calif., 2003, Vol. 54, pp. 493–529.
31. B. Fu, X. Shan, D. H. Zhang and D. C. Clary, Chem. Soc. Rev., 2017, 46, 7625–7649.
32. Y. Wang and J. M. Bowman, Chem. Phys. Lett., 2010, 491, 1–10.
33. V. Babin, G. R. Medders and F. Paesani, J. Chem. Theory Comput., 2014, 10, 1599–1607.
34. T. T. Nguyen, E. Székely, G. Imbalzano, J. Behler, G. Csányi, M. Ceriotti, A. W. Götz and F. Paesani, J. Chem. Phys., 2018, 148, 241725.
35. J. O. Richardson and S. C. Althorpe, J. Chem. Phys., 2011, 134, 054109.
36. J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso, G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate and S. C. Althorpe, Science, 2016, 351, 1310–1313.
37. J. M. Bowman, X. Wang and Z. Homayoon, J. Mol. Spectrosc., 2015, 311, 2–11.
38. C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, The MIT Press, Cambridge, Massachusetts, 2006.
39. C. Hättig, D. P. Tew and A. Köhn, J. Chem. Phys., 2010, 132, 231102.
40. O.-P. Koistinen, E. Maras, A. Vehtari and H. Jónsson, Nanosyst.: Phys., Chem., Math., 2016, 7, 925.
41. O.-P. Koistinen, F. B. Dagbjartsdóttir, V. Ásgeirsson, A. Vehtari and H. Jónsson, J. Chem. Phys., 2017, 147, 152720.
42. R. Fletcher, Practical Methods of Optimization, John Wiley and Sons, Chichester, 2nd edn, 1987.
43. D. P. Tew, W. Klopper, R. A. Bachorz and C. Hättig, in Ab Initio Theory for Accurate Spectroscopic Constants and Molecular Properties, ed. M. Quack and F. Merkt, Wiley, 2011.
44. C. Hättig, W. Klopper, A. Köhn and D. P. Tew, Chem. Rev., 2012, 112, 4–74.
45. D. P. Tew and W. Klopper, J. Chem. Phys., 2005, 123, 074101.
46. S. Ten-no, Chem. Phys. Lett., 2004, 398, 56–61.
47. E. F. Valeev, Chem. Phys. Lett., 2004, 395, 190–195.
48. D. P. Tew, W. Klopper and C. Hättig, Chem. Phys. Lett., 2008, 452, 326–332.
49. D. P. Tew and W. Klopper, J. Chem. Phys., 2006, 125, 094302.
50. D. Bokhan, S. Bernadotte and S. Ten-no, J. Chem. Phys., 2009, 131, 084105.
51. A. Köhn and D. P. Tew, J. Chem. Phys., 2010, 133, 174117.
52. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al., MOLPRO, version 2012.1, a package of ab initio programs, 2012, http://www.molpro.net.
53. T. B. Adler, G. Knizia and H.-J. Werner, J. Chem. Phys., 2007, 127, 221106.
54. D. P. Tew, J. Chem. Phys., 2016, 145, 074103.
55. B. Kolb, P. Marshall, B. Zhao, B. Jiang and H. Guo, J. Phys. Chem. A, 2017, 121, 2552–2557.
56. J. Cui and R. V. Krems, J. Phys. B: At. Mol. Phys., 2016, 49, 224001.
57. J. P. Alborzpour, D. P. Tew and S. Habershon, J. Chem. Phys., 2016, 145, 174112.
58. B. J. Braams and J. M. Bowman, Int. Rev. Phys. Chem., 2009, 28, 577–606.
59. A. P. Bartók, R. Kondor and G. Csányi, Phys. Rev. B, 2013, 87, 184115.
60. G. Montavon, K. Hansen, S. Fazli, M. Rupp, F. Biegler, A. Ziehe, A. Tkatchenko, A. von Lilienfeld and K.-R. Müller, Learning Invariant Representations of Molecules for Atomization Energy Prediction, Curran Associates, Inc., 2012, pp. 440–448.
61. F. A. Faber, A. S. Christensen, B. Huang and O. A. von Lilienfeld, J. Chem. Phys., 2018, 148, 241717.
62. S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt and K.-R. Müller, Sci. Adv., 2017, 3, e1603015.
63. P. Pulay and G. Fogarasi, J. Chem. Phys., 1992, 96, 2856–2860.
64. J. Baker, A. Kessi and B. Delley, J. Chem. Phys., 1996, 105, 192–212.
65. T. Wu, H.-J. Werner and U. Manthe, Science, 2004, 306, 2227–2229.
66. R. Welsch and U. Manthe, J. Chem. Phys., 2012, 137, 244106.
67. Y. V. Suleimanov, R. Collepardo-Guevara and D. E. Manolopoulos, J. Chem. Phys., 2011, 134, 044131.
68. K. Karandashev and J. Vaníček, J. Chem. Phys., 2015, 143, 194104.
69. A. Chakraborty, Y. Zhao, H. Lin and D. G. Truhlar, J. Chem. Phys., 2006, 124, 044315.
70. Y. Suleimanov, J. Allen and W. Green, Comput. Phys. Commun., 2013, 184, 833–840.
71. H. F. von Horsten, S. T. Banks and D. C. Clary, J. Chem. Phys., 2011, 135, 094311.
72. Q. Meng and J. Chen, J. Chem. Phys., 2017, 146, 024108.
73. D. L. Baulch, C. T. Bowman, C. J. Cobos, R. A. Cox, T. Just, J. A. Kerr, M. J. Pilling, D. Stocker, J. Troe, W. Tsang, R. W. Walker and J. Warnatz, J. Phys. Chem. Ref. Data, 2005, 34, 757–1397.
74. R. Sivaramakrishnan, J. V. Michael and B. Ruscic, Int. J. Chem. Kinet., 2012, 44, 194–205.
75. P. R. Spackman, D. Jayatilaka and A. Karton, J. Chem. Phys., 2016, 145, 104101.
76. J. O. Richardson, Phys. Chem. Chem. Phys., 2017, 19, 966–970.
77. J. O. Richardson, Faraday Discuss., 2016, 195, 49–67.
78. W. Wang and Y. Zhao, Phys. Chem. Chem. Phys., 2011, 13, 19362–19370.
79. G. N. Simm and M. Reiher, J. Chem. Theory Comput., 2017, 13, 6108–6119.
80. J. Proppe, T. Husch, G. N. Simm and M. Reiher, Faraday Discuss., 2016, 195, 497–520.
81. C. Riplinger and F. Neese, J. Chem. Phys., 2013, 138, 034106.
82. I. S. Ufimtsev and T. J. Martinez, Comput. Sci. Eng., 2008, 10, 26–34.
83. A. M. Cooper, P. P. Hallmen and J. Kästner, J. Chem. Phys., 2018, 148, 094106.

---