Calculation of reaction rate constants from a classical Wigner model based on a Feynman path integral open polymer

The Open Polymer Classical Wigner (OPCW) method [J. Chem. Phys. 152, 094111 (2020)] is applied to a flux-Heaviside trace for calculating reaction rate constants. The obtained expression for the reaction rate constant is tested on a symmetric Eckart potential. The OPCW method is shown to converge toward the exact classical Wigner result. The OPCW method shows some promise for rate constant calculations and suggestions for future tests are made. GRAPHICAL ABSTRACT


Introduction
Reaction rates are of central importance in chemistry and methods to calculate reaction rate constants are therefore of interest. In many cases classical mechanics is sufficient in order to calculate rate constants, but quantum mechanics can be essential for instance for light atoms, such as hydrogen, and in cold environments.
Methods to obtain reaction rate constants from a Wigner distribution include the classical Wigner method [1,2], with alternatives such as Feynman-Kleinert Linearised Path Integral [3] (FK-LPI), Local Gaussian Approximation Linearised Semi-Classical Initial Value Representation [4] (LGA-LSC-IVR), and to some extent the transition state theory with quantum partition functions of Wigner [5] and Eyring [6] and the quantum transition state theory of Pollak and Liao [7]. Other CONTACT G. Nyman nyman@chem.gu.se notable approximate quantum mechanical methods to calculate reaction rate constants are e.g. the instanton reaction rate theory [8][9][10], the statistical adiabatic channel model [11] and Ring Polymer Molecular Dynamics (RPMD) [12]. This article presents reaction rate constant calculations using a recent implementation of the classical Wigner method by the same authors [13]. Section 2 is an introduction to a general formulation of chemical reaction rates, then section 3 combines this theory with the Open Polymer Classical Wigner (OPCW) method. Section 4 notes some important corrections to consider for running the calculations. Sections 5, 6, and 7 respectively give the details of the computations that have been run, show the results of these computations, and contain the conclusions that have been drawn.

Reaction rate constants
Miller and coworkers published three quantum mechanical traces that can be used to acquire bimolecular reaction rate constants [14,15]. These traces are where h(x) is the Heaviside function,x is the position operator, β = (k B T) −1 , where k B is Boltzmann's constant and T is the absolute temperature,Ĥ is the Hamiltonian operator of the system, i is the imaginary unit, t is time, is the reduced Planck constant, andF is the probability density flux operator. These functions can be used to obtain reaction rate constants, k r , through the relations where Q R is the canonical partition function per unit volume of the reactants. It can be noted that the traces presented here have a symmetrised Boltzmann operator [16], i.e. e − β 2Ĥ on both sides of the first operator in each trace. The relations 4 -6 are just as valid for the real part of the trace with an asymmetrically placed Boltzmann operator [14]. Important to notice is that these relations are defined for potentials where the reactants start infinitely far from each other and the products end up infinitely far from each other. This requires the potential to be unbound. This restriction applies to the potential employed in this article. In the computations presented in this article the trace used for the calculation of rate constants is C fs (t).

Reaction rate constant calculation with the open polymer classical wigner method
In this section expressions for calculating reaction rate constants with the OPCW method will be derived. These derivations follow the ones in the original OPCW article [13] quite closely, but use the symmetrised Boltzmann operator and will be specialised to the flux-side trace, C fs (t) in Equation 2.
Initially C fs (t) is written as the integral over position, x 1 , of a matrix element.
The Boltzmann operators are each divided into N 2 factors of e − β NĤ , where N must be an even number. N−1 identity operators are inserted between these operators.
This is an imaginary time, −i β, Feynman path integral [17], which can be expressed with Wigner transforms instead of matrix elements, In the limit N → ∞ the Wigner transform of the Boltzmann operator simplifies to the classical Boltzmann factor. Also, the Boltzmann operators can be separated from the flux and Heaviside operators in this limit. We obtain where y j = x j +x (j mod N)+1 2 and H(y j , p j ) is the classical Hamiltonian. If V(y j ) is independent of momentum, all momenta except p N 2 and p N can be integrated out analytically as done in the original derivation of the OPCW method [13]. This gives where m is the mass. For the flux operatorF = 1 2m (δ(x)p +pδ(x)), the Wigner transform is where δ(x) is the Dirac delta function. Using the momentum integrals of Appendix A in our previous article [13] a quantity F (y N can be defined as This can be put into Equation 11, giving We apply the classical Wigner approximation, where x(y N , p N , t) is the position at time t of a classical trajectory starting from y N and p N . Then, by assuming finite N, we obtain which is the expression corresponding to the y 1 -version of the OPCW method for the flux-Heaviside trace. Due to the symmetrisation of the Boltzmann operator, the expression contains y N 2 rather than y 1 and below it is simply referred to as the y-version of the OPCW method. The trace will be called C fs,y (t). As shown in Appendix 1, we can also formulate an alternative version of the flux-Heaviside trace where the delta function now takes x N 2 as argument. This results in: which below is called the x-version of the OPCW method for the flux-Heaviside trace.

Correction factors for evaluation of traces with monte carlo
In this work the position integrals in Equations 16 and 17 are evaluated with Metropolis Monte Carlo [18]. In Monte Carlo integration there is a weight function and an integrand. Here, the weight function used is where y y,x is either y N 2 or x N 2 +1 . The integrand is thus for the y-version of the OPCW method and for the x-version of the OPCW method. This means that the quantity actually calculated is C fs (t) which is similar to a canonical partition function, but the integrand includes a delta function and lacks a factor of Z δ is similar to Z δ , but it has e − i p N (x 1 −x N ) in the integrand, and is thus a partition function with a delta function in it. The quantities obtained from the Monte Carlo integration are thus C fs (t) where the additional correction factor Z δ has to be computed separately.

Computational details
In this section we first define the Eckart potential. We then give some details on how the traces in Equations 16 and 17 were evaluated and how the time propagation was done. Finally we consider some numerically exact expressions against which our OPCW results can be compared.

Potential and system parameters
The symmetric Eckart barrier is given by [19] where ω is the absolute value of the angular frequency at the top of the barrier, i.e. at the transition state.

Traces
The integrals in equations 16 and 17 were evaluated through the Metropolis Monte Carlo method [18]. The maximum stepsize for the Monte Carlo was chosen and updated according to the procedure in the original OPCW work [13]. A Maxwell-Boltzmann distribution of the inverse temperature β/N was used to sample the momentum, p N . For the sampling, the ran2 pseudorandom number generator of Press et al. [20] was used. A molecular dynamics trajectory was run for each 10 th Monte Carlo step using the velocity Verlet algorithm [21,22], with a timestep of 0.05 ω −1 and a total time of 10 ω −1 .
The standard deviations of the results where calculated with the block average method, as explained in e.g. [23,24], using a minimum block size of 10 6 Monte Carlo steps.
To obtain the correction factor Z δ for the Eckart potential the numerical matrix multiplication scheme [25] was used to calculate the matrix elements of the Boltzmann operator. As an alternative, we show how to obtain Z δ by Monte Carlo in Appendix 2.

Expressions for comparison
The classical limit of the flux-Heaviside trace is easily shown to be The corresponding quantum mechanical result can be calculated from an analytic transmission coefficient   [ 19,26], where E is energy. The trace is then obtained by numerically integrating In this work the integral is performed by Romberg integration [27], using an integration interval from 0 to 100 ω and a convergence criterion of 10 −10 . The reactants are represented as a free particle, which has the partition function per unit of length This partition function is the same in the classical and quantum mechanical cases. To obtain exact classical Wigner flux-Heaviside traces for the Eckart potential the numerical matrix multiplication scheme [25] was used to calculate the matrix elements of the Boltzmann operator. An RPMD comparison was calculated using the formulation of Craig and Manolopoulos [28].

Results and discussion
In Figures 1-6 C fs (t) for the Eckart barrier at β ω = 1, 3, and 6 are shown for the x and y versions of OPCW. In Figure 7 we further show the x version result of OPCW for β ω = 12. At this low temperature, the convergence of the OPCW result with respect to the number of Monte Carlo steps becomes challenging and we chose a lower value of N to alleviate this. The rate constants and tunnelling factors, = k r,QM /k r,CM , are presented in Tables 1 and 2.
At β ω = 1 the quantum mechanical rate constant is only ∼ 6% bigger than the classical one. For the three lower temperatures the differences are significantly bigger, i.e. the quantum mechanical rate constant is ∼ 50% bigger for β ω = 3, a factor ∼ 5 bigger for β ω = 6 and finally a factor of ∼ 2000 bigger at the lowest temperature. The classical Wigner method gives rate constants that are intermediate between classical mechanics and quantum mechanics, but for the three lower temperatures they represent a significant improvement over the classical result.
The OPCW results converge to the exact classical Wigner result as the number of beads is increased. From  Figures 1-6 it can be seen that the y-version of the OPCW method needs more beads than the x-version. The xversion of the OPCW method has excellent convergence with respect to the number of beads. C fs,y (t) converge from above and C fs,x (t) converge from below. Both the x and y versions of OPCW give results that are close to the exact classical Wigner result.
The RPMD and classical Wigner methods give similar results for all temperatures. Results from a single, symmetric, one-dimensional, potential is, however, not enough to draw general conclusions.

Conclusion and outlook
In this work, the recently developed Open Polymer Classical Wigner (OPCW) method [13] is extended to the flux-Heaviside trace of Miller [15] and applied to a symmetric Eckart potential, for four temperatures. For all but the lowest temperature, the OPCW method can be converged to the exact classical Wigner result by increasing the number of beads in the path integral polymer. For the lowest temperature, a strict convergence of the OPCW result was computationally too demanding, but a result some thirty percent from the exact Classical Wigner model was obtained, see Figure 7. This should be compared to the classical prediction that is about two thousand times smaller than exact quantum mechanics.
Particularly the x-version of the OPCW method shows promise for further tests.
In the future, it would be of interest to apply the OPCW model to reaction rate constant calculations of multidimensional systems, such as the double well potential of Topaler and Makri [29] coupled to the harmonic bath of Caldeira and Leggett [30]. For such a system it could be beneficial to modify the OPCW calculations by making the bath classical, as illustrated in the preceding article [13]. There a ten-dimensional problem was considered. A more elaborate approach would be to keep an intermediate number of beads for the bath degrees of freedom, while still retaining the full number for the reaction coordinate. Such an approach will be the subject of future research. It would also be of interest to try the method on asymmetric potentials. Such potentials are challenging not just for the OPCW method but for all classical Wigner implementations, since the positioning of the dividing surface is not obvious. However, Liu and Miller have shown that this can be done successfully [4].
For systems where the reactants are bound it could be of interest to use the Heaviside-Heaviside trace of Miller et al. [14] instead of the flux-Heaviside trace.

Acknowledgments
Financial support from the Swedish Research Council (Vetenskapsrådet), diary numbers 2016-03275 and 2020-05293, is acknowledged. The computations in this work were performed at Chalmers Centre for Computational Science and Engineering (C3SE), a partner center of the Swedish National Infrastructure for Computing (SNIC).

Disclosure statement
No potential conflict of interest was reported by the author(s). To obtain an expression corresponding to the x 1 -version of OPCW we start from equation 8. There thex-dependent parts ofF have to operate to the right in the matrix element before rewriting as Wigner transforms. To be able to operate to the right, the flux operator has to be rearranged into a form where allx-dependence is to the right of allp-dependence. Using that p = −i d dx this can be achieved aŝ where the square brackets denote a commutator. Using this new form of the flux operator and operatingx to the right we obtain Rewriting as Wigner transforms, taking the limit where the Wigner transform of the Boltzmann operator is the classical Boltzmann factor, and assuming that V(y j ) is independent of momentum, a new way to write Equation 11 can be obtained: Knowing that p W x, p = p, p N 2 can be integrated out as and we get . This can be used to eliminate the derivative of the delta function through Up to this point the new derivation is equivalent to what was done for the y-version up to Equation 14. However, once this expression is subjected to the classical Wigner approximation and by assuming finite N, this is not the case anymore, as we obtain is the expression corresponding to the x 1 -version of OPCW for the flux-Heaviside trace. Due to the symmetrisation of the trace x N 2 +1 rather than x 1 is used. This trace is called C fs,x (t) in the main text.

Appendix 2
Here we show how to calculate the Z δ factor using a Monte-Carlo implementation of the effective frequency theory of Feynman and Kleinert (FK). For an introduction to the FK theory, we refer the reader to Hagen Kleinerts book [31]. The theory develops approximations to path integrals involving the Boltzmann operator. The factor Z δ is given by Equation 22: In multi-dimensional problems, it is necessary to evaluate equation A10 by Monte Carlo methods. In the FK formulation, equation A10 takes the form (see page 477 in [31])