# Assessing Parameters for Ring Polymer Molecular Dynamics Simulations at Low
Temperatures: DH+H Chemical Reaction

Ivan S. Novikov<EMAIL_ADDRESS>Yury V. Suleimanov  Alexander V.
Shapeev  Skolkovo Institute of Science and Technology, Skolkovo Innovation
Center, Nobel St. 3, Moscow 143026, Russia  Computation-based Science and
Technology Research Center, Cyprus Institute, 20 Kavafi Street, Nicosia 2121,
Cyprus

###### Abstract

Ring polymer molecular dynamics (RPMD) is an accurate method for calculating
thermal chemical reaction rates. It has recently been discovered that low-
temperature calculations are strongly affected by the simulation parameters.
Here, for the thermally activated reaction DH + H $\to$ D + H2, we calculate
the RPMD rate constants at $T=50$, 100, and 300 K and demonstrate that for
$T\geq 100$ K the standard input parameters yield accurate results, but at low
temperatures (e.g., 50 K) one must increase the asymptotic distance and force
constant, and decrease the umbrella integration step.

###### keywords:

ring polymer molecular dynamics; chemical reaction rate; quantum mechanical
effects

††journal: Chemical Physics Letters

##  1 Introduction

Accurate determination of rate constants for gas-phase bimolecular reactions
is of great importance in modeling chemical kinetics in various atmospheric,
astrochemical, and combustion systems. Chemical reactions that take place at
low temperatures are the most challenging case because quantum mechanical
effects such as zero-point energy (ZPE), tunneling, and resonance effects
become critically important at such temperatures. Recently, ring polymer
molecular dynamics (RPMD) has been proposed as an efficient tool for
calculating thermal chemical reaction rates. The method is based on the
classical isomorphism between the quantum statistical mechanics of the
physical system and the classical statistical mechanics of a fictitious ring
polymer including a number of classical copies (beads) of the original system
connected by harmonic springs [13]. RPMD is the full-dimensional approximate
quantum mechanical method for computing Kubo-transformed real-time correlation
functions for many dynamical properties, in particular, this method computes
flux-side correlation functions [33] which are used in the formalism for
calculating chemical reaction rates. The path integral nature of the RPMD
method allows capturing quantum effects, such as ZPE effects [34, 21] and
tunneling [35] that play a key role while calculating thermal rate
coefficients at low temperatures. It has been shown that the RPMD technique
[14, 15, 16, 22] provides systematically accurate approach for calculating
rate coefficients in a wide temperature range. The RPMD method has been
implemented in the RPMDrate code [39] and tested for numerous prototype
chemical reactions [34, 35, 8, 21, 31, 37, 1, 18, 38, 36, 24, 25, 26, 5].

Despite the instantaneous success of the RPMD approach, the calculation of
reaction rate constants is computationally demanding at low temperatures due
to the specifics of the technique (e.g., a large number of beads is necessary
to capture the quantum mechanical effects, long propagation times are needed
for the convergence of the rate constants, etc.). The present paper is focused
on finding the optimal set of parameters for using the RPMDrate code
effectively to predict reaction rate constants at low temperatures down to 50
K. For this aim, we selected the three-atom chemical reaction DH + H $\to$ D +
H2 which despite being simple exhibits a number of unexpected effects, such as
geometric phase effects [23, 42] and a narrow Feshbach resonance at energies
below the reaction barrier [43]. This reaction also plays an important role in
astrochemistry. The cooling of HD molecule is linked to the gravitational
collapse and the fragmentation of clouds [17, 30, 32]. In the paper [40] it
was reported that almost all deuterium can be converted to HD molecule and
subsequent cooling may lead to the formation of primordial low mass stars and
brown dwarfs.

Note that the system D+H2 has one more channel, namely, the exchange DH + H
$\to$ DH + H [3]. However, it wasn’t considered in the present study which is
focused on achieving the sufficient convergence of the RPMD thermal chemical
rate constants for the title channel D + H2.

##  2 Ring Polymer Molecular Dynamics

The ring polymer Hamiltonian for a system of $N$ atoms is written in atomic
cartesian coordinates as

$\begin{array}[]{c}\displaystyle H({\bf p},{\bf
q})=\sum_{i=1}^{N}\sum_{j=1}^{n}\left(\frac{{{\bf
p}_{i}^{(j)}}^{2}}{2m_{i}}+\frac{1}{2}m_{i}\omega^{2}\left|{\bf
q}_{i}^{(j)}-{\bf q}_{i}^{(j-1)}\right|^{2}\right)+\\\
\displaystyle\sum_{j=1}^{n}V({\bf q}_{1}^{(j)},\ldots,{\bf
q}_{N}^{(j)}),\end{array}$ (1)

where $n$ is the number of ring polymer beads, $m_{i}$ is the mass of the
$i$-th atom, ${\bf q}_{i}^{(j)}$ and ${\bf p}_{i}^{(j)}$, $j=1,\ldots,n$ are
the position and the momentum vectors of the $i$-th atom of the system,
correspondingly, ${\bf q}_{i}^{(0)}\equiv{\bf q}_{i}^{(n)}$, $i=1,\ldots,N$.
The force constant of the harmonic springs is $\omega={\beta\hbar}/{n}$, where
$\beta=1/k_{B}T$, and $T$ is the temperature of the system.

Ring polymer molecular dynamics is a classical molecular dynamics in an
extended ($n$-bead imaginary time path integral) phase space. The real-time
correlation function formalism is used for the RPMD rate coefficient
calculation. The expression for the ring polymer flux-side correlation
function [33] is given by the $t\to\infty$ limit ($t\to t_{\rm plateau}$ in
practice, where $t_{\rm plateau}$ is a “plateau” time [12] when all relevant
to the chemical event correlations decay)

$c_{\rm fs}^{(n)}(t)=\dfrac{1}{(2\pi\hbar)^{fn}}\int d^{fn}{\bf p}_{0}\int
d^{fn}{\bf q}_{0}e^{-(\beta/n)H({\bf p}_{0},{\bf q}_{0})}\delta[s({\bf
q}_{0})]\nu_{s}({\bf p}_{0},{\bf q}_{0})h[s({\bf q}_{t})],$ (2)

where $f=3N$, $s({\bf q})$ is the reaction coordinate at ${\bf q}$,
$\nu_{s}({\bf p}_{0},{\bf q}_{0})$ is the initial velocity along the reaction
coordinate

$\nu_{s}({\bf p},{\bf q})=\sum_{i=1}^{N}\sum_{j=1}^{n}\dfrac{\partial s({\bf
q})}{\partial{\bf q}_{i}^{(j)}}\dfrac{{\bf p}_{i}^{(j)}}{m_{i}},$ (3)

and $h[s({\bf q}_{t})]$ is a Heaviside step function, which counts ring
polymer trajectories $({\bf p}_{t},{\bf q}_{t})$ that are on the product side
of the dividing surface at time $t$. The RPMD rate coefficient is

$k^{(n)}(T)=\dfrac{c_{\rm fs}^{(n)}(t\to\infty)}{Q_{r}^{(n)}(T)},$ (4)

where $Q_{r}^{(n)}(T)$ is the $n$-bead path integral approximation to the
quantum mechanical reactant partition function [15, 16]. The partition
function can be quite difficult to compute accurately [7]. To avoid the direct
computation of the partition function we apply the technique implemented in
the RPMDrate code [39] and briefly described below.

The calculation of RPMD rate coefficient starts with introducing two dividing
surfaces. The first dividing surface is located in the asymptotic reactant
valley

$s_{0}(\bar{\bf q})=R_{\infty}-|\bar{\bf R}|,$ (5)

where $\bar{\bf q}=(\bar{\bf q}_{1},\ldots,\bar{\bf q}_{N})$, $\bar{\bf
q}_{i}=\dfrac{1}{n}\sum\limits_{j=1}^{n}{\bf q}_{i}^{(j)}$, $\bar{\bf R}$ is
the centroid of the vector that connects the centres of mass of the reactants,
$R_{\infty}$ is an asymptotic distance large enough to make interaction
between the reactants negligible.

The second dividing surface, $s_{1}(\bar{\bf q})$, is situated near the
transition state region and is defined in terms of the bond-breaking and bond-
forming distances. In this paper, we consider the reaction with one channel
and three atoms A, B, and C: AB + C $\to$ A + BC, i.e., one bond breaks
between the atoms A and B and one bond forms between the atoms B and C.
Therefore, the second dividing surface is determined as

$s_{1}(\bar{\bf q})=\left(|\bar{\bf q}_{\rm AB}|-q_{\rm
AB}^{\ddagger}\right)-\left(|\bar{\bf q}_{\rm BC}|-q_{\rm
BC}^{\ddagger}\right),$ (6)

where $\bar{\bf q}_{ij}$ is the vector that connects the centroids of $i$-th
and $j$-th atoms, $q_{ij}^{\ddagger}$ is the interatomic distance at the
transition state.

The next step is the introduction of the reaction coordinate $\xi$ that
connects two dividing surfaces

$\xi(\bar{\bf q})=\dfrac{s_{0}(\bar{\bf q})}{s_{0}(\bar{\bf q})-s_{1}(\bar{\bf
q})}.$ (7)

As it was mentioned above, the direct computation of the RPMD rate coefficient
with (4) can be quite difficult. Therefore, the Bennett-Chandler factorization
[4, 12] is used for the calculation:

$k_{\rm RPMD}(T)=k_{\operatorname{QTST}}(T;\xi^{\ddagger})\kappa(t\to t_{\rm
plateau};\xi^{\ddagger}).$ (8)

The first term, $k_{\operatorname{QTST}}(T;\xi^{\ddagger})$, is the quantum
transition state theory rate coefficient (its centroid-density version [20,
19, 41], cd-QTST, when the dividing surfaces are defined using the centroid
variables, as in the present case) evaluated near the transition state
$\xi^{\ddagger}$. It is defined in terms of the potential mean force (PMF), or
free energy, $W(\xi)$, along the reaction coordinate

$k_{\operatorname{QTST}}(T;\xi^{\ddagger})=4\pi
R_{\infty}^{2}\left(\frac{1}{2\pi\beta\mu_{R}}\right)^{1/2}e^{-\beta[W(\xi^{\ddagger})-W(0)]},$
(9)

where $\mu_{R}$ is the reduced mass of the reactants. The PMF difference
$W(\xi^{\ddagger})-W(0)$ is calculated using umbrella integration [28, 29, 27]

$W(\xi^{\ddagger})-W(0)=\int\limits_{0}^{\xi^{\ddagger}}\sum\limits_{i=1}^{N_{\rm
windows}}\left[\dfrac{n_{i}\mathcal{N}({\bar{\xi}}_{i},\sigma_{i}^{2})}{\sum\limits_{i=1}^{N_{\rm
windows}}n_{j}\mathcal{N}({\bar{\xi}}_{j},\sigma_{j}^{2})}\left(\dfrac{1}{\beta}\dfrac{\xi-{\bar{\xi}}_{i}}{\sigma_{i}^{2}}-k_{i}(\xi_{i}-{\bar{\xi}}_{i})\right)\right]d\xi,$
(10)

where $N_{\rm windows}$ is the number of biasing windows, $n_{i}$ is the total
number of steps sampled for the $i$-th window,
$\mathcal{N}({\bar{\xi}}_{i},\sigma_{i}^{2})$ is the normal distribution of
the reaction coordinate $\xi$ with the mean value (or, center of the biasing
window) ${\bar{\xi}}_{i}$ and the variance $\sigma_{i}^{2}$, and $k_{i}$ is
the force constant which defines the strength of the bias.

The second term in (8), $\kappa(t\to t_{\rm plateau};\xi^{\ddagger})$, is the
ring polymer recrossing factor. This factor plays a role of a dynamical
correction to cd-QTST and ensures that the resulting RPMD rate coefficient
$k_{\rm RPMD}(T)$ will be independent of the choice of the dividing surface
[14]. It is expressed as the ratio between two real-time flux-side correlation
functions

$\kappa(t\to t_{\rm plateau};\xi^{\ddagger})=\dfrac{c_{\rm fs}^{(n)}(t\to
t_{\rm plateau};\xi^{\ddagger})}{c_{\rm fs}^{(n)}(t\to
0_{+};\xi^{\ddagger})}.$ (11)

Thus, after calculating the cd-QTST rate coefficient and the recrossing
factor, we can calculate the RPMD rate coefficient.

##  3 Steps of the RPMDrate code calculation and parameters to be optimized

The RPMDrate code works in five sequential steps:

  * 1.

Generating the initial configurations for umbrella integration;

  * 2.

Umbrella integration along the reaction coordinate;

  * 3.

Calculation of the cd-QTST rate coefficient;

  * 4.

Computation of the recrossing factor;

  * 5.

Calculation of the resulting RPMD thermal rate coefficient.

We describe all the RPMDrate stages below and the crucial code parameters that
significantly affect the resulting RPMD rate constant and, thus, should be
optimized.

We start from the introduction of the umbrella sampling interval with the left
end $\xi_{\rm min}$ and the right end $\xi_{\rm max}$ and divide it into a
series of windows $\xi_{i}$, $i=1,\ldots,N_{\rm windows}$. At the first
RPMDrate step we generate classical configurations ($n=1$) for each window at
$T=300$ K. These configurations are used in subsequent umbrella integration
RPMD simulations.

At the second RPMDrate step, we equilibrate the system for each $i$-th window
(during the $t_{\rm equilibration}$ time) and accumulate ${\bar{\xi}}_{i}$ and
$\sigma_{i}$ along the RPMD thermostatted trajectory run for $t_{\rm
sampling}$. This process is repeated $N_{\rm trajectory}$ times for
accumulating enough steps of sampling and is the main part of umbrella
integration [28, 29, 27]. For the reaction investigated here, the trajectory
is thermostatted using one of the generalized Langevin equation (GLE)
thermostats [9, 10, 11]. As the result of the second RPMDrate stage, we obtain
the potential of mean force, or, free energy, $W(\xi)$, and its maximum,
$W(\xi^{\ddagger})$, which depends on temperature $T$. The cd-QTST rate
coefficient, calculated at the next RPMDrate stage, exponentially depends on
the difference between the maximum $\xi=\xi^{\ddagger}$ of free energy and its
value at $\xi=0$ (see (9), (10)) and, thus, the difference
$W(\xi^{\ddagger})-W(0)$ is the crucial value for estimation of the cd-QTST
rate coefficient. As it was demonstrated in the previous papers [39, 38, 24,
25, 26], the accuracy of the mentioned difference calculation depends on the
number of beads $n$, especially at low temperatures. Moreover, the choice of
the asymptotic distance $R_{\infty}$ and, therefore, the umbrella integration
step $d\xi$ and the force constant $k_{i}$ also affects the resulting cd-QTST
rate coefficient [5]. Thus, the parameters $n$, $R_{\infty}$, $d\xi$, and
$k_{i}$ are the most important ones to be optimized during the first three
RPMDrate steps.

Next, at the fourth RPMDrate step, we calculate the recrossing factor. For
this calculation, we use a constrained RPMD simulation in the presence of a
thermostat (parent trajectory) to generate a series of independent
configurations at the initial time moment ${\bf q}_{0}$ with centroids on the
transition state dividing surfaces $s_{\xi^{\ddagger}}(\bar{\bf
q})=\xi^{\ddagger}s_{1}(\bar{\bf q})+(1-\xi^{\ddagger})s_{0}(\bar{\bf q})=0$.
The RATTLE algorithm [2] is introduced into the time integration to constrain
the centroid $\bar{{\bf q}}_{0}$ to the dividing surfaces
$s_{\xi^{\ddagger}}(\bar{\bf q})=0$. For each of these constrained
configurations obtained from the parent trajectory, a number of momentum
vectors ${\bf p}_{0}$ is randomly sampled from the Maxwell distribution and
the resulting recrossing (child) trajectories are evolved forward in the
propagation (“plateau”) time $t_{\rm plateau}$ without the thermostat or the
dividing surface constraint. One of the reasons of inaccurate determining the
recrossing factor may be a short “plateau” time, because the decay of the
recrossing factor can occur rather slowly, and it is important to guarantee
that the recrossing factor is a constant during a long time interval. Thus,
the time $t_{\rm plateau}$ is one of the most important parameters to be
optimized at the fourth RPMDrate step.

Finally, at the fifth step we calculate the resulting RPMD thermal rate
coefficient. Note that the strategy is also suitable for the classical
dynamics (with $n=1$ beads) and at the classical high-temperature limit
($T\to\infty$).

The RPMDrate code was successfully applied to many chemical reactions. The
parameters used for calculation of different chemical reaction rates are shown
in Table 1 [34, 35, 8, 21, 31, 37, 1, 18, 38, 36, 24, 25, 26, 5]. Thermally
activated reactions were mainly considered at the temperatures $T\geq 100$ K.
For most of these reactions the standard number of beads, umbrella integration
step, and force constant were chosen: $n=128$, $d\xi=0.01$, $k_{i}=2.72$
($T/$K) eV. Thus, the above umbrella integration step and the force constant
could be chosen for calculations with any asymptotic distance $R_{\infty}$ at
the temperatures $T=100$ and 300 K. As opposed to the thermally activated
reactions, barrirless reactions were also considered at low temperatures. In
the paper Ref. [5] it was shown that while increasing the asymptotic distance
one should also increase the force constant for low temperatures. Moreover, if
the temperature of the reaction under investigation is decreased, one should
increase the number of beads for obtaining the correct rate constants.
Finally, the longer “plateau” time is needed for barrierless reactions as
compared to thermally activated ones due to slower decay of the oscillations
in the ring polymer recrossing factor.

Table 1:  Input RPMDrate code parameters for thermally activated and barrierless reactions. The force constant $k_{i}$ is reported in ($T/$K) eV, the “plateau” time $t_{\rm plateau}$ is reported in ps, and the temperature $T$ is reported in K.  Thermally activated reaction | $R_{\infty}$ | $d\xi$ | $k_{i}$ |  $n$($T$) | $t_{\rm plateau}$
---|---|---|---|---|---
Mu+H${}_{2}\to$ MuH+H [34] | $30a_{0}$ | 0.01 | 2.72 | 512(200-300) | 0.06
D+HMu $\to$ DMu+H [35] | $30a_{0}$ | 0.01 | 2.72 | 128(150-300) | 0.05
OH+H${}_{2}\to$ H+H2O [8] | $11a_{0}$ | 0.01 | 2.72 | 128(150-300) | 0.05
O(3P)+CX4 $\to$ OX+CX3 (X=H,D) [21] | $15a_{0}$ | 0.01 | 2.72 | 128(200-300) | 0.1
Cl+CX${}_{4}\to$ XCl+CX3 (X=H,D) [31] | $30a_{0}$ | 0.01 | 2.72 | 64(300) | 0.1
H+CH${}_{4}\to$ H2+CH3 [37] | $30a_{0}$ | 0.01 | 2.72 | 128(200-300) | 0.1
OH+CX${}_{4}\to$ HXO+CX3 (X=H,D) [1] | $15a_{0}$ | 0.01 | 2.72 | 128(200-300) | 0.05
X+C2H6 $\to$ HX + C2H5 (X=H,Cl,F) [18] | $15a_{0}$ | 0.01 | 2.72 | 128(200-300) | 0.05 (H)
|  |  |  |  | 0.1 (Cl)
|  |  |  |  | 0.2 (F)
Barrierless reaction | $R_{\infty}$ | $d\xi$ | $k_{i}$ |  $n$($T$) | $t_{\rm plateau}$
X(1D)+H${}_{2}\to$ XH+H (X=C,S) [38] | $20a_{0}$ | 0.01 | 2.72 | 256(50); 128(300) | 2
C+CH${}^{+}\to$ C${}_{2}^{+}$+H [36] | $15a_{0}$ | 0.01 | 2.72 | 128(20-300) | 16
C(1D)+H${}_{2}\to$ CH+H [24] | $20a_{0}$ | 0.01 | 2.72 | 256(50); 128(300) | 2
C(1D)+D${}_{2}\to$ CD+D [25] | $15a_{0}$ | 0.01 | 2.72 | 128(50); 48(300) | 2
O(1D)+X${}_{2}\to$ OX+X (X=H,D) [26] | $15a_{0}$ | 0.01 | 2.72 | 296(50); 128(300) | 2
D++H${}_{2}\to$ HD+H+ [5] | 12-30 Å | 0.01 | 2.72-20.4 | 160(20,50); 128(75,100) | 3-8

In this paper we investigate how the parameters $n$, $R_{\infty}$, $d\xi$,
$k_{i}$, and $t_{\rm plateau}$ affect the resulting RPMD rate coefficient at
different temperatures for the thermally activated reaction DH+H $\to$ D+H2.

##  4 Results and discussion

We used the BKMP2 [6] potential energy surface for all the computations
described here. First, we optimize the parameters $R_{\infty}$, $d\xi$,
$k_{i}$ at the fixed temperatures ($T=50,100$, and 300 K) and $n=64$ (we
consider this number of beads as it is typically enough for approximate
estimation the cd-QTST rate coefficient). The optimized parameters at these
temperatures are summarized in Table 2. As for the reactions in Table 1, at
the temperatures $T\geq 100$ K the standard parameters $d\xi=0.01$ and
$k_{i}=2.72$ ($T/$K) eV are suitable for computations with $10a_{0}\leq
R_{\infty}\leq 30a_{0}$, the resulting $k_{\rm QTST}$ coefficients are similar
for all temperatures of interest. However, at low temperature $T=50$ K, the
force constant increases and the umbrella integration step decreases while the
asymptotic distance is increased. Thus, at low temperatures, more attention
should be paid to the choice of $R_{\infty}$, $d\xi$, and $k_{i}$. We
demonstrate further results for the standard parameters used in the past
($R_{\infty}=15a_{0}$, $d\xi=0.01$, $k_{i}=2.72$ ($T/$K) eV) at $T=50,100$,
and 300 K, and for the modified parameters ($R_{\infty}=30a_{0}$,
$d\xi=0.005$, $k_{i}=13.6$ ($T/$K) eV) at $T=50$ K. We note, that
$R_{\infty}\leq 30a_{0}$ is sufficient for obtaining accurate enough $k_{\rm
QTST}$ at the temperature $T=50$ K, therefore, we provide the resulting RPMD
rate coefficients only for $R_{\infty}\leq 30a_{0}$.

Table 2:  Asymptotic distance $R_{\infty}$, umbrella integration step $d\xi$, and force constant $k_{i}$ optimized at $T=50,100$, and 300 K and $n=64$. The force constants are reported in ($T/$K) eV, the cd-QTST rate coefficients are reported in cm3 s-1. The standard umbrella integration step $d\xi=0.01$ and the force constant $k_{i}=2.72$ ($T/$K) eV could be chosen for calculation with any asymptotic distance at the temperatures $T\geq 100$ K. At the low temperature, $T=50$ K, the parameter $k_{i}$ increases and the step $d\xi$ decreases while increasing the distance $R_{\infty}$ for obtaining the similar cd-QTST rate coefficients. The resulting cd-QTST rate coefficient oscillates near $1.44\times 10^{-26}$ cm3 s-1, thus, the asymptotic distance $R_{\infty}=30a_{0}$ is sufficient for obtaining accurate enough rate coefficient at $T=50$ K. $T=50$ K |  $T=100$ K, $d\xi=0.01$, $k_{i}=2.72$ |  $T=300$ K, $d\xi=0.01$, $k_{i}=2.72$
---|---|---
$R_{\infty}$ | $d\xi$ | $k_{i}$ | $k_{\rm QTST}$ | $R_{\infty}$ | $k_{\rm QTST}$ | $R_{\infty}$ | $k_{\rm QTST}$
$15a_{0}$ | 0.01 | 2.72 | $1.34\times 10^{-26}$ | $15a_{0}$ | $1.19\times 10^{-23}$ | $15a_{0}$ | $5.06\times 10^{-17}$
$20a_{0}$ | 0.01 | 10.88 | $1.44\times 10^{-26}$ | $20a_{0}$ | $1.24\times 10^{-23}$ | $20a_{0}$ | $5.05\times 10^{-17}$
$30a_{0}$ | 0.005 | 13.60 | $1.46\times 10^{-26}$ | $30a_{0}$ | $1.22\times 10^{-23}$ | $30a_{0}$ | $5.09\times 10^{-17}$
$40a_{0}$ | 0.005 | 29.92 | $1.42\times 10^{-26}$ |  |  |  |

Next, we investigate the convergence of free energy profiles used to calculate
the cd-QTST rate coefficient with respect to the number of beads. The
convergence is illustrated in Fig. 1. From the plots we conclude: while
decreasing the temperatures one should increase the number of beads for the
proper estimation of the difference $W(\xi^{\ddagger})-W(0)$, and, therefore,
the cd-QTST rate coefficient. Here we choose $n=128$ at $T=300$ K, $n=256$ at
$T=100$ K, and $n=512$ at $T=50$ K. We also observe a conventional trend: free
energy barrier increases as the temperature increases.

Figure 1:  Free energy profiles for different number of beads, standard
parameters $R_{\infty}=15a_{0}$, $d\xi=0.01$, and $k_{i}=2.72$ ($T/$K) eV
($T=50,100$, and 300 K), and modified parameters $R_{\infty}=30a_{0}$,
$d\xi=0.005$, and $k_{i}=13.6$ ($T/$K) eV ($T=50$ K). The energy barrier
increases while increasing the temperature. The convergence of the energy
profiles with respect to the number of beads is observed for $n=128$ at
$T=300$ K, $n=256$ at $T=100$ K, and $n=512$ at $T=50$ K. The converged cd-
QTST rate coefficients are similar for two various sets of parameters at
$T=50$ K: $k_{\rm QTST}=4.16\times 10^{-27}$ cm3 s-1 for $R_{\infty}=15a_{0}$,
$k_{\rm QTST}=4.52\times 10^{-27}$ cm3 s-1 for $R_{\infty}=30a_{0}$ ($n=512$).

Finally, we study the convergence of the ring polymer recrossing factors with
respect to the “plateau” time at different temperatures. The recrossing
factors in the investigated temperature range are shown in Fig. 2. The longer
“plateau” time is needed for obtaining the correct ring polymer recrossing
factor if the temperature is decreased: $t_{\rm plateau}\approx 0.15$ ps at
$T=300$ K, $t_{\rm plateau}\approx 0.35$ ps at $T=100$ K, and $t_{\rm
plateau}\approx 1$ ps at $T=50$ K.

Figure 2:  Ring polymer recrossing factors for different different “plateau”
times, standard parameters $R_{\infty}=15a_{0}$, $d\xi=0.01$, and $k_{i}=2.72$
($T/$K) eV ($T=50,100$, and 300 K), and modified parameters
$R_{\infty}=30a_{0}$, $d\xi=0.005$, and $k_{i}=13.6$ ($T/$K) eV ($T=50$ K).
When the temperature of the reaction is decreased, one should increase the
“plateau” time for obtaining the correct recrossing factor. The resulting ring
polymer recrossing factor is increased if the temperature is also increased.

The resulting sets of input parameters for all the RPMDrate steps, cd-QTST
rate coefficients, recrossing factors, and RPMD rate coefficients at
$T=50,100$, and 300 K are summarized in Table 3. As it was expected for the
thermally activated reaction, the RPMD rate coefficient increases as the
temperature increases. We also note that the resulting $k_{\rm RPMD}$ rates
are similar at $T=50$ K for both the standard parameters $R_{\infty}$, $d\xi$,
$k_{i}$ and the modified ones. However, the relative deviation between them is
$\approx 10\%$. Therefore, if a very high accuracy of calculating the RPMD
rate coefficient is required one should increase the asymptotic distance at
low temperatures.

Table 3:  Input RPMDrate parameters, resulting recrossing factors and rate coefficients at $T=50,100$, and 300 K. $T$, K | $n$ | $R_{\infty}$ | $d\xi$ |  $k_{i}$, ($T/$K) eV |  $t_{\rm plateau}$, ps |  $k_{\rm QTST}$, cm3 s-1 | $\kappa$ |  $k_{\rm RPMD}$, cm3 s-1
---|---|---|---|---|---|---|---|---
300 | 128 | $15a_{0}$ | 0.01 | 2.72 | 0.15 | $5.18\times 10^{-17}$ | 0.751 | $3.89\times 10^{-17}$
300 | 128 | $20a_{0}$ | 0.01 | 2.72 | 0.15 | $5.19\times 10^{-17}$ | 0.751 | $3.90\times 10^{-17}$
300 | 128 | $30a_{0}$ | 0.01 | 2.72 | 0.15 | $5.35\times 10^{-17}$ | 0.749 | $4.01\times 10^{-17}$
100 | 256 | $15a_{0}$ | 0.01 | 2.72 | 0.35 | $1.13\times 10^{-23}$ | 0.495 | $5.59\times 10^{-24}$
100 | 256 | $20a_{0}$ | 0.01 | 2.72 | 0.35 | $1.12\times 10^{-23}$ | 0.487 | $5.45\times 10^{-24}$
100 | 256 | $30a_{0}$ | 0.01 | 2.72 | 0.35 | $1.13\times 10^{-23}$ | 0.493 | $5.57\times 10^{-24}$
50 | 512 | $15a_{0}$ | 0.01 | 2.72 | 1 | $4.16\times 10^{-27}$ | 0.271 | $1.13\times 10^{-27}$
50 | 512 | $20a_{0}$ | 0.01 | 10.88 | 1 | $4.35\times 10^{-27}$ | 0.265 | $1.15\times 10^{-27}$
50 | 512 | $30a_{0}$ | 0.005 | 13.60 | 1 | $4.52\times 10^{-27}$ | 0.269 | $1.22\times 10^{-27}$

##  5 Conclusions

In this paper, on the example of the thermally activated reaction DH + H $\to$
D + H2, we investigated the convergence of the ring polymer molecular dynamics
(RPMD) rate constants with respect to several RPMDrate code input parameters
at different temperatures (300-50 K). We demonstrated that for the
temperatures $T=100$ and 300 K the standard umbrella integration step
$d\xi=0.01$ and the standard force constant $k_{i}=2.72$ ($T/$K) eV, typically
used in the RPMDrate code, could be utilized in a wide asymptotic distance
range $10a_{0}\leq R_{\infty}\leq 30a_{0}$, the resulting RPMD rate constants
are close to each other in all the range of asymptotic distances. As opposed
to the temperatures $T=100$ and 300 K, at the low temperature $T=50$ K, the
umbrella integration step should be decreased and the force constant should be
increased (e.g., the modified input RPMDrate code parameters for
$R_{\infty}=30a_{0}$ are $d\xi=0.005$ and $k_{i}=13.6$ ($T/$K) eV) for
obtaining the correct RPMD rate constants if the asymptotic distance is
increased. The relative deviation between the resulting RPMD rates obtained
for $R_{\infty}=15a_{0}$ and $R_{\infty}=30a_{0}$ is $\approx 10\%$ at $T=50$
K whereas the same deviation is close to zero for $T\geq 100$ K. Thus, if the
RPMD rate constant should be computed very accurately then it is necessary to
carefully choose the RPMDrate code input parameters, e.g., to increase the
asymptotic distance as it was reported for the title reaction. Due to
astrochemical importance of the title reaction, lower temperatures (20 K and
below) will be the goal of our further study.

## CRediT authorship contribution statement

Ivan S. Novikov: Investigation, Visualization, Writing - Original draft, Data
Curation. Yury V. Suleimanov: Conceptualization, Methodology, Software,
Writing - Review & Editing. Alexander V. Shapeev: Conceptualization,
Resources, Writing - Review & Editing, Project administration, Funding
acquisition, Supervision.

## Declaration of Competing Interest

None.

## Acknowledgements

The work of the authors was supported by the Russian Foundation for Basic
Research (grant number 20-03-00833). Y.V.S. also acknowledges the European
Regional Development Fund and the Republic of Cyprus through the Research
Promotion Foundation (Projects: INFRASTRUCTURE/1216/0070 and Cy-Tera
NEAY$\Pi$O$\Delta$OMH/$\Sigma$TPATH/0308/31).

## References

  * [1] Joshua W Allen, William H Green, Yongle Li, Hua Guo, and Yury V Suleimanov.  Communication: Full dimensional quantum rate coefficients and kinetic isotope effects from ring polymer molecular dynamics for a seven-atom reaction OH+CH4 $\rightarrow$ CH3+H2O, 2013.
  * [2] Hans C Andersen.  Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations.  Journal of Computational Physics, 52(1):24–34, 1983.
  * [3] FJ Aoiz, M Brouard, CJ Eyles, JF Castillo, and V Sáez Rábanos.  Cumulative reaction probabilities: A comparison between quasiclassical and quantum mechanical results.  The Journal of chemical physics, 125(14):144105, 2006.
  * [4] Charles H Bennett.  Molecular dynamics and transition state theory: the simulation of infrequent events.  ACS Publications, 1977.
  * [5] Somnath Bhowmick, Duncan Bossion, Yohann Scribano, and Yury V Suleimanov.  The low temperature D++H2 $\rightarrow$ HD+H+ reaction rate coefficient: A ring polymer molecular dynamics and quasi-classical trajectory study.  Physical Chemistry Chemical Physics, 20(41):26752–26763, 2018.
  * [6] Arnold I Boothroyd, William J Keogh, Peter G Martin, and Michael R Peterson.  A refined H3 potential energy surface.  The Journal of chemical physics, 104(18):7139–7152, 1996.
  * [7] Joel M Bowman, Dunyou Wang, Xinchuan Huang, Fermin Huarte-Larrañaga, and Uwe Manthe.  The importance of an accurate ch 4 vibrational partition function in full dimensionality calculations of the H+CH4 $\rightarrow$ H2+CH3 reaction.  The Journal of Chemical Physics, 114(21):9683–9684, 2001.
  * [8] JF Castillo and YV Suleimanov.  A ring polymer molecular dynamics study of the OH+H2(D2) reaction.  Phys. Chem. Chem. Phys., 19(43):29170–29176, 2017.
  * [9] Michele Ceriotti, Giovanni Bussi, and Michele Parrinello.  Langevin equation with colored noise for constant-temperature molecular dynamics simulations.  Physical review letters, 102(2):020601, 2009.
  * [10] Michele Ceriotti, Giovanni Bussi, and Michele Parrinello.  Nuclear quantum effects in solids using a colored-noise thermostat.  Physical review letters, 103(3):030603, 2009.
  * [11] Michele Ceriotti, Giovanni Bussi, and Michele Parrinello.  Colored-noise thermostats à la carte.  Journal of Chemical Theory and Computation, 6(4):1170–1180, 2010\.
  * [12] David Chandler.  Statistical mechanics of isomerization dynamics in liquids and the transition state approximation.  J. Chem. Phys., 68(6):2959–2970, 1978.
  * [13] David Chandler and Peter G Wolynes.  Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids.  J. Chem. Phys., 74(7):4078–4095, 1981.
  * [14] Ian R Craig and David E Manolopoulos.  Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics.  J. Chem. Phys., 121(8):3368–3373, 2004.
  * [15] Ian R Craig and David E Manolopoulos.  Chemical reaction rates from ring polymer molecular dynamics.  The Journal of chemical physics, 122(8):084106, 2005.
  * [16] Ian R Craig and David E Manolopoulos.  A refined ring polymer molecular dynamics theory of chemical reaction rates.  J. Chem. Phys., 123(3):034102, 2005.
  * [17] A Dalgarno and EL Wright.  Infrared emissivities of H2 and HD.  The Astrophysical Journal, 174:L49, 1972.
  * [18] Joaquin Espinosa-Garcia, Moises Garcia-Chamorro, Jose Carlos Corchado, Somnath Bhowmick, and Yury Suleimanov.  Vtst and rpmd kinetics study of the nine-body X+C2H6 (X=H,Cl,F) reactions based on analytical potential energy surfaces.  Physical Chemistry Chemical Physics, 2020.
  * [19] M J Gillan.  Quantum-classical crossover of the transition rate in the damped double well.  J. Phys. C, 20:3621–3641, August 1987.
  * [20] M. J. Gillan.  Quantum simulation of hydrogen in metals.  Phys. Rev. Lett., 58(6):563–566, February 1987.
  * [21] Eloisa Gonzalez-Lavado, Jose C Corchado, Yury V Suleimanov, William H Green, and Joaquin Espinosa-Garcia.  Theoretical kinetics study of the O(3P)+CH4/CD4 hydrogen abstraction reaction: the role of anharmonicity, recrossing effects, and quantum mechanical tunneling.  The Journal of Physical Chemistry A, 118(18):3243–3252, 2014.
  * [22] Scott Habershon, David E. Manolopoulos, Thomas E. Markland, and Thomas F. Miller III.  Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space.  Annu. Rev. Phys. Chem., 64(1):387–413, 2013.
  * [23] Jisha Hazra, Brian K Kendrick, and Naduvalath Balakrishnan.  Geometric phase effects in ultracold hydrogen exchange reaction.  Journal of Physics B: Atomic, Molecular and Optical Physics, 49(19):194004, 2016.
  * [24] Kevin M. Hickson, Jean-Christophe Loison, Hua Guo, and Yury V. Suleimanov.  Ring-polymer molecular dynamics for the prediction of low-temperature rates: An investigation of the C(1D)+H2 reaction.  J. Phys. Chem. Lett., 6(21):4194–4199, 2015.
  * [25] Kevin M Hickson and Yury V Suleimanov.  An experimental and theoretical investigation of the C(1D) + D2 reaction.  Physical Chemistry Chemical Physics, 19(1):480–486, 2017.
  * [26] Kevin M Hickson and Yury V Suleimanov.  Low-temperature experimental and theoretical rate constants for the O (1D)+ H2 reaction.  The Journal of Physical Chemistry A, 121(9):1916–1923, 2017.
  * [27] Johannes Kästner.  Umbrella integration in two or more reaction coordinates.  J. Chem. Phys., 131(3):034109, 2009.
  * [28] Johannes Kästner and Walter Thiel.  Bridging the gap between thermodynamic integration and umbrella sampling provides a novel analysis method:“umbrella integration”.  J. Chem. Phys., 123(14):144104, 2005.
  * [29] Johannes Kästner and Walter Thiel.  Analysis of the statistical error in umbrella sampling simulations by umbrella integration.  J. Chem. Phys., 124(23):234106, 2006.
  * [30] Stephen Lepp, PC Stancil, and A Dalgarno.  Atomic and molecular processes in the early universe.  Journal of Physics B: Atomic, Molecular and Optical Physics, 35(10):R57, 2002.
  * [31] Yongle Li, Yury V Suleimanov, William H Green, and Hua Guo.  Quantum rate coefficients and kinetic isotope effect for the reaction Cl+CH4 $\rightarrow$ HCl+CH3 from ring polymer molecular dynamics.  The Journal of Physical Chemistry A, 118(11):1989–1996, 2014.
  * [32] Ian D McGreer and Greg L Bryan.  The impact of HD cooling on the formation of the first stars.  The Astrophysical Journal, 685(1):8, 2008.
  * [33] William H Miller, Steven D Schwartz, and John W Tromp.  Quantum mechanical rate constants for bimolecular reactions.  The Journal of chemical physics, 79(10):4889–4898, 1983.
  * [34] Ricardo Pérez de Tudela, F. J. Aoiz, Yury V. Suleimanov, and David E. Manolopoulos.  Chemical reaction rates from ring polymer molecular dynamics: Zero point energy conservation in Mu + H2 $\rightarrow$ MuH + H.  J. Phys. Chem. Let., 3(4):493–497, 2012.
  * [35] Ricardo Pérez de Tudela, Yury V. Suleimanov, Jeremy O. Richardson, Vicente Sáez Rábanos, William H. Green, and F. J. Aoiz.  Stress test for quantum dynamics approximations: Deep tunneling in the muonium exchange reaction D+HMu $\rightarrow$ DMu+H.  J. Phys. Chem. Lett., 5(23):4219–4224, 2014.
  * [36] Sergio Rampino and Yury V Suleimanov.  Thermal rate coefficients for the astrochemical process C+CH+ $\rightarrow$ C${}_{2}^{+}$+H by ring polymer molecular dynamics.  The Journal of Physical Chemistry A, 120(50):9887–9893, 2016.
  * [37] Yury V Suleimanov, Rosana Collepardo-Guevara, and David E Manolopoulos.  Bimolecular reaction rates from ring polymer molecular dynamics: Application to H+CH4 $\to$ H2+CH3.  J. Chem. Phys., 134(4):044131, 2011.
  * [38] Yury V Suleimanov, Wendi J Kong, Hua Guo, and William H Green.  Ring-polymer molecular dynamics: Rate coefficient calculations for energetically symmetric (near thermoneutral) insertion reactions (X+H2) $\rightarrow$ HX+H (X=C(1D),S(1D)).  The Journal of chemical physics, 141(24):244103, 2014.
  * [39] Yu.V. Suleimanov, J.W. Allen, and W.H. Green.  RPMDrate: Bimolecular chemical reaction rates from ring polymer molecular dynamics.  Comp. Phys. Comm., 184(3):833–840, 2013.
  * [40] Hideya Uehara and Shu-ichiro Inutsuka.  Does deuterium enable the formation of primordial brown dwarfs?  The Astrophysical Journal Letters, 531(2):L91, 2000.
  * [41] Gregory A. Voth, David Chandler, and William H. Miller.  Rigorous formulation of quantum transition state theory and its dynamical corrections.  J. Chem. Phys., 91:7749, 1989.
  * [42] Daofu Yuan, Yin Huang, Wentao Chen, Hailin Zhao, Shengrui Yu, Chang Luo, Yuxin Tan, Siwen Wang, Xingan Wang, Zhigang Sun, et al.  Observation of the geometric phase effect in the H+ HD $\rightarrow$ H2\+ D reaction below the conical intersection.  Nature communications, 11(1):1–7, 2020.
  * [43] Boyi Zhou, Benhui Yang, Naduvalath Balakrishnan, Brian K Kendrick, and Phillip C Stancil.  Prediction of a feshbach resonance in the below-the-barrier reactive scattering of vibrational excited HD with H.  The Journal of Physical Chemistry Letters, 2020.