Accelerating quantum instanton calculations of the kinetic isotope effects

Path integral implementation of the quantum instanton approximation currently belongs among the most accurate methods for computing quantum rate constants and kinetic isotope effects, but its use has been limited due to the rather high computational cost. Here we demonstrate that the efficiency of quantum instanton calculations of the kinetic isotope effects can be increased by orders of magnitude by combining two approaches: The convergence to the quantum limit is accelerated by employing high-order path integral factorizations of the Boltzmann operator, while the statistical convergence is improved by implementing virial estimators for relevant quantities. After deriving several new virial estimators for the high-order factorization and evaluating the resulting increase in efficiency, using $\mathrm{\cdot H_{\alpha}+H_{\beta}H_{\gamma}\rightarrow H_{\alpha}H_{\beta}+\cdot H_{\gamma }}$ reaction as an example, we apply the proposed method to obtain several kinetic isotope effects on $\mathrm{CH_{4}+\cdot H\rightleftharpoons\cdot CH_{3}+H_{2}}$ forward and backward reactions.


I. INTRODUCTION
Accurate evaluation of the rate constant, i.e., the prefactor of the rate law of elementary chemical reactions, remains one of the central goals of chemical kinetics because this constant reflects the mechanism of the reaction as well as other properties of the potential energy surface (PES) on which the reaction occurs. Another quantity that is frequently used for studying reaction mechanisms, and, in particular, detecting nuclear quantum effects on reaction rates, is the kinetic isotope effect (KIE). The KIE is defined as the ratio of rate constants for two isotopologs, i.e., molecules that only differ in isotope composition. These effects, which include, e.g., tunneling, corner-cutting, and zero-point energy effect, tend to play an important role in hydrogen transfer reactions with a high activation barrier. Although they are most important at low temperatures, nuclear quantum effects sometimes manifest themselves even at physiological temperatures, a fact uncovered by studying KIE's on some enzymatic reactions. 1 Several approaches are currently used for calculating rate constants and KIE's in situations where quantum effects are not negligible. One approach consists in adding a tunneling correction to transition state theory, 2 others use an approximation for the propagator by treating it semiclassically 3 or by treating only one or two degrees of freedom quantum mechanically. 4 Another promising method is the ring polymer molecular dynamics 5 (RPMD), which can partially capture both quantum effects and classical recrossing. Finally, there are various quantum generalizations of the transition state theory. These so-called quantum transition state theories 6,7 include the quantum instanton (QI) approximation to the rate constant, 8 i.e., the method whose efficiency we attempt to increase in the present paper. The a) Electronic mail: jiri.vanicek@epfl.ch QI approximation is motivated by the semiclassical instanton theory [9][10][11][12] and, as the name suggests, takes into account only the zero-time properties of the reactive flux-flux correlation function; however, in contrast to the semiclassical instanton, the QI approximation treats the Boltzmann operator exactly quantum mechanically. This improvement makes QI quite accurate as verified in many previous applications of the method. [13][14][15][16][17][18] QI theory expresses the reaction rate in terms of imaginary-time correlation functions, which, in turn, can be evaluated by path integral (PI) Monte Carlo (MC) methods. 19 As for KIE's, the problem can be simplified further by using thermodynamic integration. 20 The resulting method, however, has a drawback common to all PI methods: it operates in a configuration space of greatly increased dimensionality, leading to high computational cost. Indeed, the quantum limit is approached as the number of dimensions goes to infinity. Several approaches have been proposed to bypass the problem and this paper combines two of them to accelerate the QI calculations.
The first approach employs Boltzmann operator factorizations of higher order of accuracy. The resulting PI representations of relevant quantities exhibit faster convergence to the quantum limit, allowing to reduce the dimensionality of the calculation. [21][22][23][24][25][26] The second approach uses improved estimators with lower statistical errors, which permit shortening the MC simulation. 27,28 In this work, we combine these two strategies and, in addition, propose several new estimators. We then test the resulting method on two systems: the model ·H α + H β H γ → H α H β + ·H γ rearrangement, for which we also evaluate the resulting gain in computational efficiency, and the reaction CH 4 + ·H ·CH 3 + H 2 , a process whose KIE's were studied in detail both experimentally and theoretically, with classical transition state theory (TST), several of its corrected versions, 29,30 reduced dimensionality quantum dynamics, 31 and RPMD. 32 TABLE I. Summary of the notation used in this paper for a system of N particles in a D-dimensional Euclidean space. m i is the mass of particle i; v and w are vectors defined in the N D-dimensional configuration space, while v i and w i are their D-dimensional components corresponding to particle i; A is a Hermitian matrix defined over the configuration space, and A i j is its D × D dimensional submatrix containing only the columns corresponding to particle i and rows corresponding to particle j.

Expression Comment
∇ i Gradient with respect to coordinates of particle i v i · w i  D j=1 v i, j · w i, j Standard dot product of v i and w i in the D-dimensional Mass-weighted inner product of v and w in the system's configuration space, where s ∈ {−1, 0, 1} on the right-hand side, while on the left-hand side, a corresponding shorthand notation s ∈ {−, 0, +} is used.
(The value of s depends on whether v or w are covariant or contravariant.) ∥v∥ s  ⟨v, v⟩ s Mass-weighted norm of a configuration space vector ⟨v, A, w⟩ su Matrix product of A with v and w; the same shorthand notation is used for s and u as in ⟨v, w⟩ s The rest of the paper is organized as follows: After outlining the derivation of the QI approximation for the rate constant in Sec. II, in Sec. III, we first show how this approximation can be combined with the PI formalism and then explain in detail the two strategies to improve numerical performance of the standard PI implementation. The numerical results are presented in Sec. IV, while Sec. V concludes the paper. To facilitate the reading, our notation is summarized in Table I.

II. QUANTUM INSTANTON FORMALISM
The QI approximation for the thermal rate constant k(T) can be derived from the exact Miller-Schwartz-Tromp formula, 11 expressing the product of the rate constant with the reaction partition function Q r as the time integral of the flux-flux correlation function where is the symmetrized correlation function of operatorsÂ andB at temperature T = 1/(k B β) and time t, is the operator of flux through dividing surface (DS) γ ∈ {a, b}, r is the position vector in the N D-dimensional configuration space (N is the number of atoms, D is the number of spatial dimensions), and h is the Heaviside function [i.e., h(x) = 1 for x ≥ 0 and h(x) = 0 for x < 0]. The two DS's a and b completely separate the reactant and product regions and are defined by the equation ξ γ (r) = 0. In addition, ξ γ are chosen so that ξ γ (r) > 0 for r in the product region and ξ γ (r) < 0 in the reactant region.
The QI approximation can be derived by applying the steepest descent approximation to Eq. (1). 20,33 First, one multiplies and divides the integrand of Eq. (1) by the so-called delta-delta correlation function where ∆ γ is the normalized delta function, and ∥ · ∥ − is the norm of a covariant vector (see Table I).
Then, one assumes that C ff (t) decays sufficiently fast so that the main contribution to the integral in Eq. (1) comes from t close to zero (hence the name "quantum instanton"), and that for these short times, the ratio C ff (t)/C dd (t) remains approximately constant and given by C ff (0)/C dd (0). One can therefore evaluate the time integral in Eq. (1) with the steepest descent approximation, obtaining the QI expression for the rate constant where is a certain energy variance. For reasons that will become clear below, we keep C dd (0) in Eq. (8), even though it may seem to cancel out. The last question to be addressed is how to choose positions of the optimal DS's. From semiclassical considerations, it follows that the best choice is to require that C dd (0) be a saddle point with respect to ξ a and ξ b ; 8 if ξ γ are controlled by a set of parameters {η (γ) k }, the stationarity condition becomes ∂C dd

III. GENERAL PATH INTEGRAL IMPLEMENTATION
The QI approximation allows expressing the rate constant in terms of the reactant partition function and properties of flux-flux and delta-delta correlation functions at time t = 0. In this section, we first explain how the PI formalism allows transforming the quantum problem of finding these quantities to a classical one, applied to the so-called polymer chain, 34,35 and then describe an efficient implementation allowing a significant acceleration of calculations of the KIE's.
One of our goals is using higher-order factorizations of the Boltzmann operator in order to accelerate the convergence of the KIE's to the quantum limit. In Subsection III A, we therefore present a general derivation of the PI expression for the Boltzmann operator matrix element, valid for all Boltzmann operator factorizations used in this work, and in Subsection III B, we obtain general PI expressions for Q r and C dd (0). In Subsections III C and III D, we explain how all quantities necessary for computing the KIE within the QI approximation can be expressed in terms of thermodynamic averages over ensembles corresponding to PI expressions for Q r and C dd (0); in Subsection III E, we derive estimators allowing to calculate these averages with a lower statistical error and therefore significantly accelerating statistical convergence, which is our second main goal. Some of the more tedious derivations are deferred to Appendix A.

A. Lie-Trotter (LT), Takahashi-Imada (TI), and Suzuki-Chin (SC) factorizations of the imaginary-time path integral propagator
The coordinate matrix element of the Boltzmann operator at temperature T = 1/(k B β) can be rewritten as a matrix element of the product of P ∈ N Boltzmann operators at a higher temperature inversely proportional to the parameter ϵ β/P, We next consider three possible high-temperature factorizations of the Boltzmann operator as follows: 1. The symmetrized version of the Lie-Trotter factorization, This second-order factorization, which we will, for simplicity, denote by LT is the one most commonly used for discretizing the imaginary-time Feynman PI. 2. The TI factorization, 36 Tr is an effective one-particle potential. This fourth-order factorization significantly accelerates the convergence to the quantum limit of the PI expression for Q r . However, it only behaves as a fourth-order factorization when it is used for evaluating the trace of the Boltzmann operator. If one naively removes the Tr in Eq. (13) and applies the resulting factorization e −ϵĤ ≈ e −ϵV TI /2 e −ϵT e −ϵV TI /2 (15) to off-diagonal elements, which are required for PI representations of C dd (0) and C ff (0), one obtains only secondorder convergence, and no numerical advantage over the LT factorization. Since it will allow us to provide a single derivation of many quantities for different factorizations, we will abuse terminology and refer to Eq. (15) also as "Takahashi-Imada" factorization, keeping in mind that the original authors were aware that their splitting is of the fourth order only in the context of Eq.
are the "endpoint" and "midpoint" effective one-particle potentials. The dimensionless parameter α can assume an arbitrary value, but evidence in the literature 21,23 suggests that α = 0 gives superior results in most PI simulations, and hence, it was also the value used in our calculations. Now we use one of the three PI splittings for each of the P high-temperature factors in Eq. (11), with the caveat that for the SC factorization (only), we replace P with P/2 (so P must be even) and ϵ = β/P with ϵ = 2 β/P in Eq. (11). After inserting (P − 1) resolutions of identity in the coordinate basis in front of every kinetic factor (except the first one), we obtain where the effective potentialΦ and prefactor C are defined as In the expression forΦ, we use the notation r (P) r (b) , r (0) r (a) for the boundary points; N is the number of atoms, D is the number of spatial dimensions, m i is the mass of particle i, ∥ · ∥ + is the norm of a contravariant vector (see Table I), and V (s) eff is the effective one-particle potential, where is the coordinate representation of the commutator term in Eqs. (14), (17), and (18). In the context of discretized PI's, the integer P is often referred to as the Trotter number. The coefficient d s for the fourth-order correction of an effective one-particle potential depends on the splitting used, LT splitting, 1/24, TI splitting, α/6, SC splitting and s even, (1 − α)/12, SC splitting and s odd.
The weightsw s in the sum over effective one-particle potentials also depend on the splitting: for endpoint s (i.e., s = 0, P), these weights arew s = 1/2 for the LT and TI splittings, andw s = 1/3 for the SC splitting; for other values of s,w s = 1 for the LT and TI splittings, whereas for the SC splitting,w s = 4/3 for odd s andw s = 2/3 for even s. Expression (19) becomes exact as P goes to infinity.

B. Path integral representations of the partition and delta-delta correlation functions
From Eq. (19), it is straightforward to obtain the PI representation Q r, P of the reactant partition function Q r ; in particular, (In general, we will distinguish between a quantity A and its PI representation A P for a given value of P by adding an additional subscript P.) By  d{r (s) }, we mean integration over all r (s) , s ∈ {1, . . . , P}; ρ r ({r (s) }) can be regarded as the thermal distribution of {r (s) } of the new system; Φ is the closed-loop version ofΦ, i.e.,

Φ({r (s) })
Φ (r (P) , r (1) , . . . , r (P) ) From now on, we will always consider closed loops such that r (0) = r (P) . The difference between the new weights w s and the old weightsw s is that w P = 1 for the LT or TI splittings and w P = 2/3 for the SC splitting, for which we also require P to be even.
We can now see that the PI representation of the quantum partition function Q r is identical to the classical partition function Q r,cl of a system (called "polymer chain") where every original particle is replaced with P pseudoparticles connected by harmonic forces. Also note that for P = 1 and the LT factorization, our PI expression reduces to the expression for the classical Boltzmann distribution.
For C dd (0), we have, analogously, (we shall omit the time argument of C dd and C ff if it equals 0). Note that ρ ‡ ({r (s) }) differs from ρ r ({r (s) }) by the two delta constraints imposed on r (P/2) and r (P) .
In the rest of the section, we will show how the QI expression for the KIE can be rewritten in terms of classical thermodynamic averages over ρ r and ρ ‡ . Expressions for the corresponding estimators will be presented in a general way valid for all Boltzmann operator splittings considered in this work and as such will contain the main part common for all splittings and a part which corresponds to the fourth-order corrections and is only non-zero if a splitting other than LT is used; since this additional part depends on the gradient of the potential energy, we will denote it by adding "grad" subscript to the name of the estimator. Although it is one of the main results of this work, for clarity, the derivation of the parts associated with the fourth-order factorizations will be left for Appendix A. Before we proceed, it is necessary to point out relative costs of running MC simulations over ρ r and ρ ‡ obtained with different Boltzmann operator splittings. While the use of the LT splitting only requires potential energy for each r (s) , the SC splitting with 0 < α < 1 and the TI factorization also require the gradient of energy for each r (s) , and the SC splittings with α = 0 and α = 1 require gradients for r (s) with s odd and even, respectively.

C. Estimators for constrained quantities
Within the PI formalism, both the energy spread ∆H and the flux factor C ff /C dd can be expressed as thermodynamic averages over the ensemble whose configurations are weighted by ρ ‡ ({r (s) }). 19 In order to obtain the PI representation of C dd (0) and ∆H 2 , it is convenient to perform the Wick rotation and define a new function of a complex argument ζ. Supposing that C dd (t) is analytic, The PI representation of C dd (ζ) is with prefactor C = P and two "partial" effective potentials wherew s =w s for all s except for s = P/2, for which w P/2 =w P =w 0 . The effective potentialsΦ + andΦ − in Eq. (33) are obtained in a similar manner asΦ was obtained in Eq. (19). The difference is that instead of the matrix element of the Boltzmann operator exp(− βĤ), one considers an element of exp(− β +Ĥ /2) or exp(− β −Ĥ /2), and the exponential operators are discretized into P/2 rather than P parts. As a result, expressions (37) and (38) forΦ + andΦ − can be obtained from the one forΦ [Eq. (20)] if β is replaced with β + /2 and β − /2, respectively, and P is replaced with P/2. After differentiating expression (33) with respect to ζ, using Eq. (32) to go from C ′′ dd, P (ζ) back toC dd, P (t), and noting that After the substitution of expressions (29) and (39) for C dd, P andC dd, P into definition (9) of ∆H 2 , the estimator for ∆H 2 takes the form is used as the weight function. From now on, if a quantity A can be expressed as a classical thermodynamic average, we will denote the corresponding estimator by A est (the density function over which the averaging is performed will not be denoted explicitly since this will always be clear from the context). Explicit differentiation in Eqs. (40) and (41) leads to the so-called thermodynamic estimator, 19 The ratio C ff /C dd can be computed by the Metropolis algorithm as well. To obtain the corresponding estimator, we first note that the flux operator can be expressed aŝ CombiningF γ with the PI representation of the Boltzmann operator, one obtains 19 where f v is the so-called velocity factor, r a = r (P/2) , r b = r (P) , ⟨·, ·⟩ − is the inner product of two covariant vectors, and ⟨·, ·, ·⟩ −− the matrix product of a covariant matrix with two covariant vectors (see Table I).
Taking the ratio of PI representations (46) and (29) of C ff and C dd immediately yields the estimator for the ratio C ff /C dd , The thermodynamic estimator takes the form 19 where ⟨·, ·⟩ 0 is the inner product of a covariant and contravariant vectors (see Table I), and we defined

D. Thermodynamic integration with respect to mass
The last ingredient needed for evaluating QI rate constant (8) is the ratio C dd /Q r , which, unfortunately, cannot be calculated by the standard Metropolis algorithm. However, in the case of KIE's, one can circumvent this problem by employing the so-called thermodynamic integration with respect to mass, 20 which is easy to understand from the explicit QI expression for the KIE, where A and B are different isotopologs of otherwise the same system. The basic idea of the thermodynamic integration with respect to mass consists in computing ratios r by considering a continuous transformation 20,39-41 from A to B using a dimensionless parameter λ ∈ [0, 1] controlling atomic masses of the intermediate systems as Ratios r are rewritten in terms of their logarithmic derivatives, which are normalized quantities and, therefore, can be calculated with the Metropolis algorithm, The integrals in the exponent can be evaluated numerically with Simpson's rule or other standard methods. However, several approaches have been proposed for decreasing the exponentiated integration error in the ratio Q (A) r /Q (B) r , which can further accelerate the calculation by lowering the required number of integration points: these include rescaling of mass (which in the simplest variant involves 25,41 or introducing higher-order derivatives of Q r with respect to λ; 25 if the ratio is close to unity, it is also possible to eliminate the integration error altogether by using direct estimators for 42 In the case of d ln C dd /dλ needed in the QI, one needs to keep track of the possible change of ξ γ during the course of the integration, In Ref. 20, the authors proposed to choose {η (γ) k (λ)} that satisfy Eq. (10) at each λ integration step, making the second term in Eq. (56) exactly zero and leaving only ∂ ln C dd /∂λ to be considered. Here, we take an alternative and more numerically stable approach: By introducing new accurate estimators for ∂ ln C dd /∂η (γ) k , we can avoid having to find the optimal values of η (γ) k (λ) for all λ. Instead, we only find optimal η (γ) k (λ) at the boundary points λ ∈ {0, 1}, obtain other, not necessarily optimal, η (γ) k (λ) by linear interpolation, and evaluate both terms of Eq. (56) for each λ.
The estimators for ∂ ln Q r /∂λ and ∂ ln C dd /∂λ are where F ds is the contribution that comes from differentiating the mass-dependent normalization factor in Eq. (6), Here, ∇ i is the gradient with respect to coordinates of particle i (see Table I).
Direct evaluation of Eqs. (57) and (58) yields the thermodynamic estimators 20 Derivation of the estimator for ∂ ln C dd /∂η (γ) k involves a rather tedious algebra and is therefore presented in Appendix B; the result is where ∇ (γ) is the gradient with respect to r γ and is the term associated with the change of configuration space volume satisfying the constraint. Obtaining the thermodynamic estimator for B k(γ) is straightforward and yields

E. Virial estimators
So far we have only considered thermodynamic estimators, which are obtained via direct differentiation of the Boltzmann operator matrix elements. However, an estimator for a given quantity is not unique; it is often possible to obtain an estimator with smaller statistical error. Among such estimators are the so-called virial and centroid virial estimators, 43,44 which are motivated by the virial theorem of classical mechanics and which can be derived 22,45 most simply by applying a coordinate transformation before the differentiation.
Two of the five virial estimators used in this work, namely, the estimators for ∂ ln Q r /∂λ and ∆H 2 had been proposed previously; 27,28 the former, however, had not been used in combination with the SC factorization. To derive the centroid virial estimator for ∂ ln Q r /∂λ, let us choose an arbitrary bead u and rewrite Q r in terms of the coordinates where {m ′ i } are a set of parameters with dimensionality of mass. One then substitutes the new C and Φ resulting from the transformation of coordinates into Eq. (57), and, finally, sets {m ′ i } = {m i } and transforms back to initial coordinates. This procedure yields an improved virial estimator, which, however, depends on an arbitrary choice of bead u. After taking the arithmetic average of all P estimators corresponding to P different choices of u ∈ {1, . . . , P}, the centroid virial estimator is obtained, where is the centroid coordinate. (Explicit form of F r,grad is derived in Appendix A.) From now on, we will refer to this estimator as "virial"; originally, the name "centroid virial" was introduced to distinguish the estimator from the simple virial estimator derived in Ref. 28, which was not considered in this work since its statistical error is larger than the error of its centroid counterpart.
For ∆H 2 , one starts 27 by rewriting Eq. (33) using the coordinates whereř (s) is the reference point given by The kinetic parts ofΦ ± are rewritten in the new coordinates; for example, forΦ + , one uses the relation By substituting transformedΦ ± and C into Eqs. (40) and (41), one obtains the desired G and F terms of the virial estimator, where ⟨·, ·, ·⟩ 00 is the matrix product of a covariant matrix with two contravariant vectors (see Table I), and explicit forms of G v,grad and F grad are derived in Appendix A. Now let us turn to the derivation of the new estimators promised in the Introduction. In particular, we propose new virial estimators for ∂ ln C dd /∂λ, C ff /C dd , and ∂ ln C dd /∂η (γ) k . For ∂ ln C dd /∂λ, we use a coordinate rescaling which is similar to Eq. (69) and yields the virial estimator For C ff /C dd , we introduce new coordinates and employ the identity RewritingΦ ± in terms of {x (1) , . . . , x (P/2−1) , r (P/2) , x (P/2+1) , . . . , x (P−1) , r (P) } and inserting them into Eq. (47) lead to the virial estimator where we introduced coefficients Using the same rescaling as for f v,v , we can also derive the virial estimator for ∂ ln C dd /∂η (γ) k , where r (P/2) γ stands for r (P) if γ = a and for r (P/2) if γ = b. We would like to comment on the cost of using the estimators described in this subsection. While thermodynamic estimators require little numerical effort, their virial counterparts depend on the gradient and Hessian of the effective potential. (Note that although B k(γ) th also depends on the force, it depends only on the force acting on a single bead, and hence, its cost is negligible for large P.) It should be emphasized, however, that gradient-and Hessian-dependent parts of virial estimators can be calculated by finite difference, without the need to evaluate the gradient or Hessian explicitly. For example, ⟨w, ∇V ⟩ 0 and ⟨w, ∂ 2 V/∂r 2 , w⟩ 00 are the first and second derivatives of V in the direction of w and therefore can be evaluated by finite difference using just one and two additional evaluations of V , respectively. As a result, the effective cost is only one extra potential evaluation per bead for F r,cv , one per unconstrained bead for F ‡ cv , two per unconstrained bead for (G v + F 2 v )/2, and three for f v,v . Calculating B k(γ) th will require exactly one potential evaluation and calculating B k(γ) v will require P − 1 evaluations unless it is computed at the same time as f v,v (in this case, it would require just one extra potential evaluation, other numerical ingredients being shared with f v,v ).
It should be emphasized that it is not necessary to evaluate these estimators after each MC step due to finite correlation lengths inherent to MC simulations. This realization frequently allows one to make the additional cost of evaluating even the more expensive estimators small compared with the cost of the random walk itself.
Finally, we would like to point out that, while authors of Refs. 22 and 28 used finite differences with respect to mass and β, respectively, to calculate virial estimators of interest, we found this approach less convenient since it requires introducing two parameters (finite difference steps) that must be adjusted for each new isotopolog and for each temperature. We therefore only used finite differences with respect to coordinates in the system's configuration space, with a single finite difference step which is the same for all isotopologs and all temperatures.

IV. APPLICATIONS
In summary, to compute the KIE on a reaction, one must do as follows: 1. Estimate the Trotter number P that is sufficient to adequately describe the system. For this purpose, we made several preliminary calculations to estimate the P necessary for the lowest and highest temperature; for other temperatures, we used the empirical rule that 1/P stays approximately linear with respect to T.
2. Choose the two DS's. We chose ξ γ (r) of the form For reactions where atom X breaks its bond with atom Y and forms a bond with atom Z, we used as a reaction coordinate the difference of the "bond" lengths, i.e., where R XY is the distance between X and Y . Optimal values of η ‡ γ were found by running test simulations to find the sign of ∂ ln C dd /∂η ‡ γ at different values of (η ‡ a , η ‡ b ). 3. Run simulations at different values of λ in order to obtain the corresponding logarithmic derivatives of Q r and C dd , as well as C ff (0)/C dd (0) and ∆H for λ = 0 and λ = 1, then evaluate Eqs. (54) and (55) using, e.g., Simpson's rule. For many systems, d ln C dd /dλ and ∂ ln Q r /∂λ are quite smooth functions and nine intermediate points were sufficient to accurately evaluate the thermodynamic integrals (i.e., the discretization error (DE) of the λ integral was smaller than the already small statistical error). After this, evaluating the KIE using Eq. (52) is straightforward.
For each value of λ, one has to run two MC simulations in {r (s) }: a "constrained simulation" with two slices constrained to their respective DS and a standard ("unconstrained") simulation. Since treating exact constraints is not straightforward in MC methods, we approximated the delta constraint with a "smeared" delta function δ sm , In contrast with the approximation used in Ref. 19, the width σ of our Gaussian δ m does not depend on temperature or mass. The approximate constraint converges to the exact delta function as σ → 0. Presence of δ sm [ξ γ (r)] can be easily simulated by adding an extra constraining potential to two of the slices. For MC sampling, we employed the staging algorithm [46][47][48] with multislice moves in combination with whole-chain moves. For constrained simulations, we also made extra single-slice moves of slices P/2 and P, since these slices are more rigid than others due to the presence of the constraining potential.

A. H + H 2 rearrangement
The errors of PI MC calculations come mostly from two sources: the PI discretization error (due to P being finite) and the statistical error inherent to MC methods. (As for quantities evaluated with thermodynamic integration, there is an additional discretization error of the thermodynamic integral due to taking a finite number of λ steps.) To verify the improvements outlined in Sec. III, we studied their influence on the behavior of the two main types of

Computational details
Statistically converged simulations (paralleled over 64 trajectories, 4 × 10 7 MC steps each) were run with different values of the Trotter number (from P = 8 to 64 with step 4 and from 64 to 352 with step 16) and different Boltzmann operator factorizations. Virial estimators were evaluated only after every 25 MC steps, whereas the thermodynamic -after every step, because the additional cost was negligible. To estimate statistical errors of the results, we calculated root mean square deviations of averages over different trajectories. [Having a relatively high number (64) of uncorrelated trajectories, we could thus avoid a more tedious block-averaging procedure, 50 but we did check in several cases that the two approaches gave very similar statistical error estimates.] As for the positions of the DS's, for calculating the KIE, choosing η ‡ a = η ‡ b = 0 was quite satisfactory even at T = 200 K (in this case, C dd is stationary from symmetry considerations) for analyzing numerical behavior of ∂ ln C dd /∂λ, ∆H 2 , and C ff /C dd . For ∂ ln C dd /∂η ‡ a , however, we used η ‡ a = −0.5 and η ‡ b = 0.5 in order to make the logarithmic derivative statistically relevant.
For this particular setup, the increase of central processing unit (CPU) time associated with evaluating all virial estimators at once was about 15% for constrained and 3.5% for unconstrained simulations. The increase of CPU time associated with the use of higher-order splittings was negligible for constrained simulations; for unconstrained simulations, it was 2.5% and 5% for SC and TI splittings, respectively.

Results
Convergence of different quantities to their quantum limits as a function of the Trotter number P is shown in Fig. 1. As expected, the SC factorization allows to lower the Trotter number significantly in comparison with the LT factorization. In the case of ∂ ln Q r /∂λ, the SC splitting is slightly outperformed by the TI factorization, which has a smaller prefactor of the error term, possibly because the TI splitting leads to an expression invariant under cyclic bead permutations.
Statistical errors of different estimators are presented in Fig. 2. Note that they do not depend much on the factorization used. In contrast, the decrease of statistical errors associated with using virial estimators is remarkable for all quantities.
To compare the speedups achieved by different combinations of splittings and estimators, we estimated the relative CPU times needed to converge the quantities ∆H, C ff /C dd , r /Q ( A) r , and C (B) dd /C ( A) dd to 1% discretization and statistical errors. 51 To estimate the speedup associated with calculating the overall KIE itself with 1% statistical and discretization errors, we ran a separate set of simulations with λ = 1 in addition to those for λ = 0; the statistical and discretization errors of the KIE calculated with different combinations of estimators and factorizations were then approximated with the corresponding errors obtained if thermodynamic integration of Q r and C dd had Results correspond to the KIE ·H + H 2 / · D + D 2 at 200 K. "v" stands for "virial" and "th" for "thermodynamic." been performed using a single step trapezoidal rule (i.e., based just on the two boundary points λ = 0 and λ = 1).
Let us assume the CPU time of a simulation to be approximately proportional to P and the number of MC steps. Then, for a given combination of factorization and estimator, the cost of achieving the target discretization and statistical errors is proportional to the productPσ 2P , whereP is the value of the Trotter number that yields the target discretization error and σP is the statistical error exhibited by the estimator at this value of P. These estimates of CPU cost are then corrected by the increase in CPU time associated with using the fourth-order splittings and virial estimators. The final results are presented in Table II, which confirms that the combination of virial estimators and fourth-order splittings leads to a significant speedup of the calculation.
One may be surprised that the value of P necessary to achieve 1% convergence of C ff /C dd appears to be roughly independent of the splitting used; this is probably because the discretization errors of C dd and C ff cancel to a larger extent for the LT than the SC splitting. Taking the discretization error to be 0.5% rather than 1% makes the difference in the required value of P even more pronounced: P = 40 for the LT and P = 80 for the SC splitting.
Note that even though the values of P required to converge individual quantities are quite large (up to P = 336 for ∂ ln Q r /∂λ if LT splitting is used), the Trotter number P necessary to converge the final KIE result is significantly lower due to the cancellation of discretization errors between individual quantities and especially between the two isotopologs. However, our P value required for the   53,54 factorizations for finding different quantities associated with the RPMD expression for the reaction rate. The authors found that for dynamical properties, the TI splitting gives little improvement over the standard LT factorization, which is consistent with our explanation presented in Subsection III A; both factorizations are outperformed by the fourth-order Chin factorization, which is in agreement with the SC outperforming LT splitting in Table II. For equilibrium properties, the authors found that the efficiencies of the Chin and TI factorizations are similar, and that both fourth-order factorizations significantly outperform the standard LT splitting, again in agreement with our results and explanation.
We mentioned earlier that we had calculated virial estimators by finite difference, making the computational cost of their evaluation independent of dimensionality. To employ fourth-order splittings, however, one must know the potential gradient for all P replicas (for the TI splitting) or at least for P/2 replicas (for the SC splitting if α = 0 or α = 1). In general, if evaluating the gradient becomes too expensive compared to the potential energy itself, it may be advantageous to use the LT instead of the fourth-order splittings. For example, as shown in Table II, using the fourth-order splittings decreased the necessary P approximately four times; therefore, for this particular system, it is reasonable to use the TI factorization if the cost of evaluating the gradient is smaller than three times the cost of evaluating the potential alone. For the SC factorization, the corresponding factor is around six, since one needs only P/2 force evaluations. This upper bound for efficiency may be pushed further using the reweighting-based techniques; 21-23 this approach, however, is known to increase the statistical errors of the final result in high-dimensional systems. 55 Finally, we verified the modified methodology by comparing our result for the KIE ·H + H 2 / ·D + D 2 with those of Ref. 20, obtained both with the QI approximation and with an exact quantum method. For each temperature, we calculated ∆H and C ff /C dd by PIMC simulations with 1.28 × 10 8 steps at λ = 0 and λ = 1. Ratios of Q r and C dd were evaluated by rewriting them as in Eqs. (54) and (55), respectively, and finding the integral over λ using Simpson's rule with integration step ∆λ = 0.1. At T = 200 K, we also ran calculations with ∆λ = 0.05 to verify that the integration error of the final result is lower than the statistical error. Values of ∂ ln C dd /∂λ within the integration interval were obtained by running simulations with 6.4 × 10 7 MC steps (i.e., fewer steps than for the λ-endpoint simulations because ∂ ln C dd /∂λ and ∂ ln Q r /∂λ tend to converge faster than C ff /C dd and especially than ∆H). These conditions ensured that the total relative error of the final KIE caused by statistical noise was below 1%. We chose P in such a way that the relative error due to P being finite was less than the statistical one. At the lowest temperature T = 200 K, we chose P = 64, while for  T = 2400 K, P = 12 turned out to be appropriate; for other temperatures, we estimated the necessary P by interpolation assuming that the 1/P is a linear function of T. To verify that the chosen values of P were sufficient, we ran additional simulations at temperatures 200 K, 1000 K, and 2400 K with λ = 0 and λ = 1 with a doubled value of P. If two KIE's calculated with ∆λ = 1 at the two different values of P differed by a value that was lower than the sum of their statistical errors, the lower value of P was deemed sufficient for the calculation. The statistical errors, i.e., root mean square errors (RMSE), were estimated with the "block-averaging" method 50 in order to remove the effect of correlation length of the random walk in the Metropolis MC simulation. In Ref. 20, η ‡ γ were taken to be 0 for all temperatures and all values of λ. Even though this choice of DS positions leads to C dd being stationary, it is a local minimum rather than a saddle point. We therefore also checked the result for the case when the proper optimal DS positions are found. Since from symmetry considerations the optimal DS parameters satisfy η ‡ a = −η ‡ b , simple bisection was sufficient to calculate the values up to 0.01 a.u. The results are presented in Table III. Intermediate results of the calculations are presented separately in Table XII in Appendix C. We can see that the values obtained with η ‡ γ = 0 agree well with those of Ref. 20, validating our modifications. It can also be seen that the full DS optimization improves agreement of the QI results with the exact quantum result, making the method remarkably accurate at low temperatures.

B. CH 4 + ·H ·CH 3 + H 2 exchange
As mentioned, the KIE's on the CH 4 + ·H ·CH 3 + H 2 exchange had been studied by various numerical methods, but not by the QI approximation. We therefore decided to test the accelerated QI method on this reaction using the potential energy surface published in Ref. 56.

Computational details
We first ran a series of trial simulations to roughly determine the value of P and the number of MC steps ensuring that at the lowest temperature, the relative statistical error of the KIE is below 1% and that the discretization error with respect to P is even smaller. The target statistical error was guaranteed by running 6.4 × 10 7 step MC simulations at λ = 0 and λ = 1, and 3.2 × 10 7 simulations at other values of λ. The target discretization error was achieved with P = 80 for the LT and P = 20 for the combination of fourth-order splittings at T = 400 K; at other temperatures, P was chosen such that the ratio β/P stayed approximately constant. We chose ∆λ = 0.1 as for the case of ·H + H 2 / · D + D 2 ; to be completely sure that the thermodynamic integration error was negligible to the statistical one, we also ran calculations with ∆λ = 0.05 at T = 400 K for the equilibrium isotope effect ·CH 3 / · CD 3 and KIE ·CH 3 + D 2 / · CD 3 + D 2 , as these cases exhibited the most drastic changes of properties during thermodynamic integration.
To determine the stationary positions of the DS's ({η ‡ γ }), we ran several short (8 × 10 6 steps) simulations to find the sign of ∂ ln C dd /∂η ‡ γ at different DS positions; the saddle points were found with accuracy of 0.01 a.u. The difference between η ‡ a and η ‡ b turned out to be negligible at all temperatures considered, in accordance with what is expected at "high" temperatures. 13 The calculated values of η ‡ are presented in Table XIII in Appendix C; as expected, they are   quite close to the position of the classical transition state at η ‡ = −0.94.

Results
Next, we compared results obtained by the accelerated method (employing a combination of fourth-order splittings and virial estimators) and by the standard method (employing a combination of LT splitting and thermodynamic estima-tors). The corresponding numerical results are labeled as "accel." and "std.," respectively. For further comparison, we calculated the same KIE's also with the conventional TST [57][58][59] and TST with Wigner tunneling correction 60 (in the tables, the corresponding columns are denoted as "TST" and "TST + Wigner", respectively). In the TST framework, the expression for the rate constant takes the form  where Q ‡ and Q r are partition functions of the transition and reactant states, computed assuming separability of rotations and vibrations, harmonic approximation for vibrations, and rigid rotor approximation for rotations. Note that the usual factor exp(−E a /k B T), where E a is the activation energy, is absorbed into our definition of Q ‡ since we use the same zero of energy for both Q ‡ and Q r . This expression can be multiplied by the so-called Wigner tunneling correction to account for tunneling contribution to the reaction rate. Here, ω ‡ is the imaginary frequency corresponding to the motion along the reaction coordinate, where µ ‡ is the effective reduced mass of the movement along the reaction coordinate at the saddle point, and K ‡ is the corresponding negative force constant. Since the conventional TST expression captures the changes of zero point energy as well as of the rotational and translational partition functions due to the isotopic substitution, one may expect that the difference between the QI and conventional TST should be largely due to the difference between the extent of tunneling present in the two isotopologs. The results are presented in Tables IV-X. First of all, it can be seen that for KIE's due to mass changes not affecting the transferred atom (see Tables IV-VI), the QI values are close to those obtained by conventional TST. This can be understood qualitatively from expression (86) for Wigner tunneling correction for reaction rates. The main contribution to µ ‡ appearing in the expression for ω ‡ comes from the transferred atom; therefore, if its mass does not change, the Wigner tunneling corrections for different isotopologs will have similar values and largely cancel out in the KIE.
Second, note that, in agreement with the usual difference in magnitudes of secondary and primary isotope effects, replacing ·CH 3 with ·CD 3 leads to a much smaller rate change than does replacing H 2 with D 2 (compare Tables V and VI and IX and X). This consideration also explains why the KIE's corresponding to ·CH 3 + H 2 / ·CH 3 + D 2 and ·CD 3 + H 2 / ·CD 3 + D 2 are quite close to each other (see Tables X and IX). For some KIE's presented in Tables VIII-X, it appears that results obtained with TST or TST with Wigner tunneling correction are in better agreement with experimental values than those obtained with the QI, probably indicating that a large cancellation takes place between the errors of the TST and of the PES.
In order to estimate the influence of the used force field on the final result, we also ran calculations with the PES published in Ref. 67 for CH 4 + ·H/CH 4 + ·D. After finding the optimal DS positions (see Table XIII in Appendix C), we compared the QI values of this KIE obtained with the two PES's from Refs. 67 and 56 (see Table XI), finding that the choice of the PES affects the KIE value by as much as 10%. In contrast, comparison of the KIE's computed with the same PES, but with two different accurate quantum methodologies (RPMD and QI), results in a remarkable agreement, within the statistical error of less than 2%. Finally, note that the QI KIE is in much better agreement with experiment if computed with the PES of Ref. 56 than with the PES of Ref. 67, suggesting that the former PES, which was used for most of the calculations in this paper, was the appropriate choice.
As for the performance of the fourth-order splittings, since an analytical gradient was not available for the CH 4 + ·H system, the gradient had to be calculated numerically using finite differences. For constrained simulations, this made the force twelve times (once per each internal degree of freedom) as expensive as the potential itself, leading to a seven-fold increase in CPU time for a given P and number of MC steps when the SC splitting was used. Since the fourth-order splitting decreased the necessary P by a factor of four, the final increase in CPU time for a given discretization error and number of MC steps was 75%. For unconstrained simulations employed to find ·CH 3 / · CD 3 equilibrium isotope effect, the force was six times as expensive as the potential; since the use of the TI factorization allowed to decrease P four times, the final increase in CPU time was also 75% for a given number of MC steps and discretization error.
In summary, the KIE's were reproduced in a reasonable agreement with experiment. The differences are probably due to both the error of the potential energy surface used and the large experimental error. Note that our accelerated methodology again drastically reduced both the discretization and statistical errors of the calculations.

V. CONCLUSIONS
In conclusion, we have accelerated the methodology from Ref. 20 for computing KIE's with the QI approximation. In particular, we have combined virial estimators (several of which have been derived for the first time here) with highorder factorizations of the quantum Boltzmann operator, and shown that this combination significantly accelerates the QI calculations of the KIE's in systems with prominent quantum effects. We have also proposed and demonstrated the utility of a new method for the thermodynamic integration of the deltadelta correlation function C dd , which is a convenient alternative to the approach employed in Ref. 20

APPENDIX A: DERIVATION OF THE FOURTH-ORDER CORRECTIONS FOR DIFFERENT ESTIMATORS
When one of the fourth-order factorizations is used, V (s) eff (r (s) ) has an explicit dependence on mass and β; as a result, one needs to add appropriate "corrections" to the estimators arising from the differentiation with respect to these quantities.
For ∂ln Q r /∂λ/ β and ∂ln C dd /∂λ/ β, it follows from Eqs. (57) and (58) that the correction F r,grad is Note that when a coordinate rescaling is used to obtain an estimator (e.g., for centroid virial estimators), the correction remains the same due to the following equality: As for F th and F v , since the gradient correction to be added is Again, this correction is the same for the virial and thermodynamic variants.
Since the G factor involves the second derivatives with respect to β, the corrections will be different for G th and G v . While G th,grad is obtained easily as to find G v , one needs to take advantage of the following relations: The only terms for which the explicit β dependence plays a role are the last three. As a result, we get This expression can be rewritten as w s d s 3V grad +  (r (s) −ř (s) ), ∇V grad (r (s) )  0  . (A10)

APPENDIX B: DERIVATION OF B k (γ)
The present derivation is a slight generalization of the one found in Ref. 68. We start by transforming to mass-scaled coordinates, which will simplify the subsequent algebra due to the equality In mass-scaled coordinates, PI representation (29) of the delta-delta correlation function can be rewritten as where the second normalized delta function has been absorbed into ρ in order to simplify the following derivation. Differentiation of C dd, P with respect to the DS's parameters yields After integrating by parts with respect to x γ in the first integral, we get Equation (62) is obtained by substituting the explicit expression for ρ and transforming back to Cartesian coordinates.

APPENDIX C: ADDITIONAL NUMERICAL RESULTS
In this section, we present some additional numerical results that were moved from the main text for the sake of clarity. Figure 3 depicts the logarithmic plots of the discretization errors of various ingredients of the QI approximation as functions of the Trotter number P. The discretization error for a quantity A is defined as | A P − A ∞ |, where A ∞ was estimated by averaging A P over several highest values of P, for which the discretization error was considered negligible. The averaging was performed in order to reduce the statistical error. The plots in Fig. 3 demonstrate the faster convergence to the quantum limit achieved with higher-order factorizations: indeed, especially for the logarithmic derivative of Q r , one can see that the discretization error dependence approaches the asymptotic behavior O(P −2 ) for the LT and O(P −4 ) for the SC and TI factorizations. In addition, in all panels, it is clear for which value of P the discretization error becomes smaller than the statistical error, since for higher values of P, the smooth dependence of the discretization error on P is obscured by statistical noise.    Table III. All quantities as well as their statistical errors are in atomic units.