Accelerating quantum instanton calculations of the kinetic isotope effects

TL;DR

This paper enhances quantum instanton calculations of kinetic isotope effects by combining high-order path integral factorizations and virial estimators, significantly reducing computational costs and improving efficiency.

Contribution

The authors develop new virial estimators for high-order path integral factorizations, enabling faster and more accurate quantum instanton calculations of kinetic isotope effects.

Findings

Achieved orders-of-magnitude efficiency increase in calculations.

Successfully applied method to methane isotope reactions.

Demonstrated improved convergence to quantum limits.

Abstract

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.

Tables

Table 1. Table 1: Summary of the notation used in this paper for a system of N 𝑁 N particles in a D 𝐷 D -dimensional Euclidean space. m i subscript 𝑚 𝑖 m_{i} is the mass of particle i 𝑖 i ; 𝐯 𝐯 \mathbf{v} and 𝐰 𝐰 \mathbf{w} are vectors defined in the N ​ D 𝑁 𝐷 ND -dimensional configuration space, while v i subscript v 𝑖 \mathrm{v}_{i} and w i subscript w 𝑖 \mathrm{w}_{i} are their D 𝐷 D -dimensional components corresponding to particle i 𝑖 i ; 𝐀 𝐀 \mathbf{A} is a Hermitian matrix defined over the configuration space, and A i ​ j subscript A 𝑖 𝑗 \mathrm{A}_{ij} is its D × D 𝐷 𝐷 D\times D dimensional submatrix containing only the columns corresponding to particle i 𝑖 i and rows corresponding to particle j 𝑗 j .

Expression	Comment
${\nabla }_i$	gradient with respect to coordinates of particle $i$
$v_i\cdot w_i:={\sum }_{j=1}^D{v}_{i,j}\cdot w_{i,j}$	standard dot product of $v_i$ and $w_i$ in the $D$-dimensional Euclidean space
${⟨𝐯,𝐰⟩}_s:={\sum }_{i=1}^N{\left(m_i\right)}^s{v}_i\cdot w_i$	mass-weighted inner product of $𝐯$ and $𝐰$ in the system’s configuration space, where $s\in \{-1,0,1\}$ on the right-hand side, while on the left-hand side a corresponding shorthand notation $s\in \{-,0,+\}$ is used
${\Vert 𝐯\Vert }_s:=\sqrt{{⟨𝐯,𝐯⟩}_s}$	mass-weighted norm of a configuration space vector
${⟨𝐯,𝐀,𝐰⟩}_{su}:={\sum }_{i=1}^N{\sum }_{j=1}^N{m}_i^s{m}_j^u{v}_i\cdot A_{ij}\cdot w_j$	matrix product of $𝐀$ with $𝐯$ and $𝐰$; the same shorthand notation is used for $s$ and $u$ as in ${⟨𝐯,𝐰⟩}_s$

Table 2. Table 2: Estimated speedups of the QI calculations achieved by the use of various combinations of path integral factorizations and estimators (th = thermodynamic, v = virial). Speedup “1” (i.e., no speedup) corresponds to the standard method employing a combination of the Lie-Trotter factorization and thermodynamic estimators. Results correspond to the KIE ⋅ H + H 2 / ⋅ D + D 2 \mathrm{\cdot H+H_{2}/\cdot D+D_{2}} at 200 K.

Factorization		Lie-Trotter (LT)				Suzuki-Chin (SC)				Takahashi-Imada (TI)
Estimator		$P$	th	v		$P$	th	v		$P$	th	v
$\partial \ln Q_r/\partial \lambda$		336	$1$	$220$		$96$	$10$	$850$		64	44	1200
$\partial \ln C_{\text{dd}}/\partial \lambda$		$240$	$1$	$19$		$64$	$13$	$160$		_	_	_
$\Delta H^2$		$128$	$1$	$170$		$48$	$14$	$380$		_	_	_
$C_{\text{ff}}/C_{\text{dd}}$		$36$	$1$	$2.7$		$32$	$1.3$	$3.0$		_	_	_
$\partial \ln C_{\text{dd}}/\partial {\eta }_a^{‡}$		$80$	$1$	$1.7$		$16$	$3.6$	$3.8$		_	_	_
KIE		$128$	$1$	$34$		$40$	$12$	$97$111All quantities except for $\partial \ln Q_r/\partial \lambda$ are calculated with the SC factorization. For $\partial \ln Q_r/\partial \lambda$ TI factorization is used.		_	_	_

Table 3. Table 3: Kinetic isotope effect ⋅ H + H 2 / ⋅ D + D 2 \mathrm{\cdot H+H_{2}/\cdot D+D_{2}} at different temperatures.

$T$ (K)	optimal ${\eta }_a^{‡}=-{\eta }_b^{‡}$		QI		QM111Ref. Vaníček et al., 2005. QM denotes exact quantum-mechanical results from this reference.	% error222The error is defined as $\left({\mathrm{KIE}}_{\mathrm{QI}}-{\mathrm{KIE}}_{\mathrm{QM}}\right)/{\mathrm{KIE}}_{\mathrm{QM}}\times 100\%$ for the optimized DS case.	QI111Ref. Vaníček et al., 2005. QM denotes exact quantum-mechanical results from this reference.
	$\lambda =0$	$\lambda =1$	no DS optimization	optimized DS
200	$1.00$	$0.41$	$22.3\pm 0.2$	$22.6\pm 0.3$	$22.53$		$23.15$
250	$0.62$	$0$	$10.91\pm 0.08$	$9.92\pm 0.09$	$10.40$	$-5$	$10.98$
300	$0.01$	$0$	$7.38\pm 0.05$	$7.35\pm 0.05$	$6.97$	$5$	$7.41$
400	$0$	$0$	_	$4.87\pm 0.03$	$4.74$	$3$	$4.84$
600	$0$	$0$	_	$3.29\pm 0.02$	$3.42$	$-4$	$3.25$
1000	$0$	$0$	_	$2.23\pm 0.01$	$2.61$	$-15$	$2.22$
1500	$0$	$0$	_	$1.81\pm 0.01$	$2.27$	$-20$	$1.83$
2400	$0$	$0$	_	$1.55\pm 0.01$	_	_	$1.56$

Table 4. Table 4: Kinetic isotope effect CH 4 + ⋅ H / CH 4 + ⋅ D \mathrm{CH_{4}+\cdot H}/\mathrm{CH_{4}+\cdot D} .

$T$ (K)	TST	TST +Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$111Ref. Pu and Truhlar, 2002	RDQD222Ref. Kerkeni and Clary, 2004	RPMD333Ref. Li et al., 2013	Expt.444Ref. Kurylo, Hollinden, and Timmons, 1970
			accel.	std.
400	0.56	0.56	0.60$\pm 0.01$	0.62$\pm 0.07$	0.54	0.58	0.64		0.74555Values taken from Ref. Pu and Truhlar, 2002
500	0.66	0.66	0.70$\pm 0.01$	0.65$\pm 0.08$	0.65	0.67	1.03	0.65	0.84555Values taken from Ref. Pu and Truhlar, 2002
600	0.74	0.74	0.78$\pm 0.01$	0.7$\pm 0.1$	0.73	0.74	1.23		0.91555Values taken from Ref. Pu and Truhlar, 2002
700	0.79	0.79	0.84$\pm 0.01$	0.9$\pm 0.1$	0.78	0.79	1.33	0.80	0.97555Values taken from Ref. Pu and Truhlar, 2002

Table 5. Table 5: Kinetic isotope effect ⋅ CH 3 + D 2 / ⋅ CD 3 + D 2 \mathrm{\cdot CH_{3}+D_{2}}/\mathrm{\cdot CD_{3}+D_{2}} .

$T$ (K)	TST	TST+Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$111Ref. Pu and Truhlar, 2002	Expt.222Based on data from Refs. Shapiro and Weston, 1972; Tsang and Hampson, 1986; Kerr and Parsonage, 1976
			accel.	std.
400	0.73	0.74	0.76$\pm 0.01$	0.7$\pm 0.1$	0.75	0.74	0.59333Values taken from Ref. Truhlar et al., 1992
500	0.82	0.82	0.83$\pm 0.01$	0.9$\pm 0.1$	0.83	0.82	0.72333Values taken from Ref. Truhlar et al., 1992
600	0.87	0.88	0.88$\pm 0.01$	0.8$\pm 0.2$	0.88	0.88	0.82333Values taken from Ref. Truhlar et al., 1992
700	0.91	0.91	0.90$\pm 0.01$	0.8$\pm 0.2$	0.92	0.91	0.90333Values taken from Ref. Truhlar et al., 1992

Table 6. Table 6: Kinetic isotope effect ⋅ CH 3 + H 2 / ⋅ CD 3 + H 2 \mathrm{\cdot CH_{3}+H_{2}}/\mathrm{\cdot CD_{3}+H_{2}} .

$T$ (K)	TST	TST+Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$111Ref. Pu and Truhlar, 2002	Expt.222Ref. Shapiro and Weston, 1972
			accel.	std.
400	0.74	0.74	0.80$\pm 0.01$	0.78$\pm 0.08$	0.75	0.81	0.85333Values taken from Ref. Espinosa-García and Corchado, 1996
500	0.82	0.83	0.86$\pm 0.01$	0.9$\pm 0.1$	0.83	0.88	0.86333Values taken from Ref. Espinosa-García and Corchado, 1996
600	0.87	0.88	0.90$\pm 0.01$	1.0$\pm 0.2$	0.88	0.92	0.87333Values taken from Ref. Espinosa-García and Corchado, 1996
700	0.91	0.91	0.92$\pm 0.01$	1.0$\pm 0.2$	0.92	0.95	0.88333Values taken from Ref. Espinosa-García and Corchado, 1996

Table 7. Table 7: Kinetic isotope effect ⋅ CH 3 + HD / ⋅ CH 3 + DH \mathrm{\cdot CH_{3}+HD}/\mathrm{\cdot CH_{3}+DH} .

$T$ (K)	TST	TST+Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$111Ref. Pu and Truhlar, 2002	Expt.222Ref. Shapiro and Weston, 1972
			accel.	std.
467	1.51	1.86	2.10$\pm 0.02$	2.5$\pm 0.4$	1.50	1.83	2.1$\pm$0.5333Values taken from Ref. Espinosa-García and Corchado, 1996
531	1.48	1.76	1.84$\pm 0.02$	2.0$\pm 0.3$	1.47	1.71	1.9$\pm$0.3333Values taken from Ref. Espinosa-García and Corchado, 1996
650	1.44	1.64	1.59$\pm 0.02$	1.4$\pm 0.3$	1.43	1.56	1.2$\pm$0.3333Values taken from Ref. Espinosa-García and Corchado, 1996

Table 8. Table 8: Kinetic isotope effect ⋅ CD 3 + HD / ⋅ CD 3 + DH \mathrm{\cdot CD_{3}+HD}/\mathrm{\cdot CD_{3}+DH} .

$T$ (K)	TST	TST+Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$111Ref. Pu and Truhlar, 2002	Expt.222Ref. Shapiro and Weston, 1972
			accel.	std.
400	1.56	1.99	2.52$\pm 0.02$	2.3$\pm 0.3$	1.55	1.91	1.85333Values taken from Ref. Pu and Truhlar, 2002
500	1.50	1.80	1.95$\pm 0.02$	1.9$\pm 0.3$	1.49	1.60	1.61333Values taken from Ref. Pu and Truhlar, 2002
600	1.46	1.68	1.65$\pm 0.02$	1.3$\pm 0.3$	1.45	1.56	1.47333Values taken from Ref. Pu and Truhlar, 2002
700	1.43	1.60	1.52$\pm 0.01$	1.6$\pm 0.3$	1.42	1.49	1.38333Values taken from Ref. Pu and Truhlar, 2002

Table 9. Table 9: Kinetic isotope effect ⋅ CD 3 + H 2 / ⋅ CD 3 + D 2 \mathrm{\cdot CD_{3}+H_{2}}/\mathrm{\cdot CD_{3}+D_{2}} .

$T$ (K)	TST	TST+Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$111Ref. Pu and Truhlar, 2002	Expt.222Ref. Shapiro and Weston, 1972
			accel.	std.
400	3.45	4.39	5.60$\pm 0.04$	5.0$\pm 0.8$	3.22	4.13	3.33333Values taken from Ref. Pu and Truhlar, 2002
500	2.98	3.57	3.92$\pm 0.03$	3.9$\pm 0.5$	2.83	3.21	2.88333Values taken from Ref. Pu and Truhlar, 2002
600	2.64	3.04	3.15$\pm 0.03$	2.6$\pm 0.6$	2.54	2.73	2.61333Values taken from Ref. Pu and Truhlar, 2002
700	2.40	2.68	2.75$\pm 0.02$	2.5$\pm 0.4$	2.33	2.43	2.43333Values taken from Ref. Pu and Truhlar, 2002

Table 10. Table 10: Kinetic isotope effect ⋅ CH 3 + H 2 / ⋅ CH 3 + D 2 \mathrm{\cdot CH_{3}+H_{2}}/\mathrm{\cdot CH_{3}+D_{2}} .

$T$ (K)	TST	TST+Wigner	QI		TST111Ref. Pu and Truhlar, 2002	CVT/$\mu \mathrm{OMT}$ 111Ref. Pu and Truhlar, 2002	Expt.222Ref. Shapiro and Weston, 1972
			accel.	std.
400	3.45	4.41	5.93$\pm 0.05$	5.8$\pm 0.8$	3.22	4.57	4.8$\pm$0.4333Values taken from Ref. Espinosa-García and Corchado, 1996
500	2.97	3.58	4.09$\pm 0.04$	4.0$\pm 0.6$	2.83	3.43	3.5$\pm$0.2333Values taken from Ref. Espinosa-García and Corchado, 1996
600	2.64	3.05	3.21$\pm 0.03$	3.1$\pm 0.5$	2.54	2.86	2.8$\pm$0.2333Values taken from Ref. Espinosa-García and Corchado, 1996

Table 11. Table 11: Influence of the potential energy surface (PES) on the KIE CH 4 + ⋅ H / CH 4 + ⋅ D \mathrm{CH_{4}+\cdot H}/\mathrm{CH_{4}+\cdot D} . Comparison of the QI KIE’s calculated using the PES’s of Refs. Corchado, Bravo, and Espinosa-Garcia, 2009 and Zhang, Braams, and Bowman, 2006 . Note also the remarkable agreement between the KIE’s computed with RPMD and QI on the same PES.

$T$ (K)	PES of Ref. Corchado, Bravo, and Espinosa-Garcia,2009				PES of Ref. Zhang, Braams, and Bowman,2006	Expt.222Ref. Kurylo, Hollinden, and Timmons, 1970
	TST	TST+Wigner	RPMD111Ref. Li et al., 2013	QI	QI
400	0.52	0.52		0.54$\pm 0.01$	0.60$\pm 0.01$	0.74 333Values taken from Ref. Pu and Truhlar, 2002
500	0.63	0.63	0.65	0.64$\pm 0.01$	0.70$\pm 0.01$	0.84 333Values taken from Ref. Pu and Truhlar, 2002
600	0.71	0.71		0.73$\pm 0.01$	0.78$\pm 0.01$	0.91 333Values taken from Ref. Pu and Truhlar, 2002
700	0.77	0.77	0.80	0.79$\pm 0.01$	0.84$\pm 0.01$	0.97 333Values taken from Ref. Pu and Truhlar, 2002

Table 12. Table 12: Values of the factors entering the QI expression ( 52 ) for the KIE ⋅ H + H 2 / ⋅ D + D 2 \mathrm{\cdot H+H_{2}/\cdot D+D_{2}} with optimized dividing surface positions, displayed in Table  IV.1.2 . All quantities as well as their statistical errors are in atomic units.

$T$ (K)	$\Delta H^2\times {10}^6$		$C_{\mathrm{ff}}/C_{\mathrm{dd}}\times {10}^3$		$C_{\mathrm{dd}}$ ratio	$Q_r$ ratio
	$\lambda =0$	$\lambda =1$	$\lambda =0$	$\lambda =1$
200	$3.68\pm 0.03$	$4.27\pm 0.02$	$1.33\pm 0.01$	$2.03\pm 0.01$	$43.8\pm 0.3$	$1404\pm 1$
250	$4.87\pm 0.04$	$4.95\pm 0.03$	$2.01\pm 0.01$	$2.10\pm 0.01$	$55.8\pm 0.3$	$572.9\pm 0.3$
300	$7.40\pm 0.04$	$5.10\pm 0.03$	$2.73\pm 0.01$	$1.97\pm 0.01$	$49.5\pm 0.1$	$316.1\pm 0.2$
400	$7.92\pm 0.05$	$6.27\pm 0.04$	$2.49\pm 0.01$	$2.01\pm 0.01$	$33.96\pm 0.06$	$149.5\pm 0.1$
600	$12.0\pm 0.1$	$10.7\pm 0.1$	$2.77\pm 0.01$	$2.44\pm 0.01$	$23.11\pm 0.04$	$70.92\pm 0.03$
1000	$26.8\pm 0.1$	$24.5\pm 0.1$	$3.82\pm 0.03$	$3.55\pm 0.02$	$18.15\pm 0.02$	$39.42\pm 0.01$
1500	$54.3\pm 0.2$	$50.3\pm 0.2$	$5.29\pm 0.04$	$5.05\pm 0.02$	$16.78\pm 0.02$	$30.11\pm 0.01$
2400	$124.6\pm 0.4$	$117.1\pm 0.4$	$7.98\pm 0.05$	$7.80\pm 0.03$	$16.36\pm 0.01$	$25.51\pm 0.01$

Table 13. Table 13: Optimal positions of the dividing surfaces along the reaction coordinate [see Eq. ( 83 )] for transition states of several isotopic variants of the CH 4 + ⋅ H ⇌ ⋅ CH 3 + H 2 \mathrm{CH_{4}+\cdot H}\rightleftharpoons\mathrm{\cdot CH_{3}+H_{2}} exchange at several temperatures.

Potential energy surface of Ref. Zhang, Braams, and Bowman,2006
TS	400 K	500 K	600 K	700 K
$H_{3}C\cdots H\cdots H$	-0.91	-0.90	-0.88	-0.86
$H_{3}C\cdots H\cdots D$	-0.88	-0.87	-0.85	-0.84
$D_{3}C\cdots H\cdots D$	-0.89	-0.88	-0.86	-0.85
$D_{3}C\cdots D\cdots H$	-0.93	-0.91	-0.89	-0.86
$H_{3}C\cdots D\cdots D$	-0.89	-0.87	-0.86	-0.85
$D_{3}C\cdots D\cdots D$	-0.90	-0.89	-0.87	-0.85
$D_{3}C\cdots H\cdots H$	-0.92	-0.90	-0.89	-0.87
	467 K	531 K	650 K
$H_{3}C\cdots H\cdots D$	-0.87	-0.86	-0.84
$H_{3}C\cdots D\cdots H$	-0.90	-0.89	-0.87
Potential energy surface of Ref. Corchado, Bravo, and Espinosa-Garcia,2009
TS	400 K	500 K	600 K	700 K
$H_{3}C\cdots H\cdots H$	-1.03	-1.00	-0.97	-0.95
$H_{3}C\cdots H\cdots D$	-1.00	-0.97	-0.94	-0.92

Equations

$$k(T)Q_r={\int }_0^{\infty }{C}_{\text{ff}}(t)dt\text{,}$$

$$C_{\text{ff}}(t):=C_{{\hat{F}}_a{\hat{F}}_b}(t)\text{,}$$

$$C_{\hat{A}\hat{B}}(t):=\mathrm{Tr}\left(\hat{A}{e}^{-(\beta /2-\mathrm{i}t)\hat{H}}\hat{B}{e}^{-(\beta /2+\mathrm{i}t)\hat{H}}\right)$$

$${\hat{F}}_{\gamma }:=-\frac{\mathrm{i}}{\hslash }[h[{\xi }_{\gamma }(\hat{\mathbf{r}})],\hat{H}]$$

$$C_{\text{dd}}(t):=C_{{\hat{\Delta }}_a{\hat{\Delta }}_b}(t)\text{,}$$

$${\Delta }_{\gamma }(\mathbf{r})=\delta [{\xi }_{\gamma }(\mathbf{r})]||\nabla {\xi }_{\gamma }(\mathbf{r})|{|}_{-}$$

$${\int }_0^{\infty }{C}_{\text{ff}}(t)dt$$

$$\approx \frac{C_{\text{ff}}(0)}{C_{\text{dd}}(0)}{\int }_0^{\infty }{C}_{\text{dd}}(t)dt$$

$$\approx \frac{C_{\text{ff}}(0)}{C_{\text{dd}}(0)}\frac{\hslash \sqrt{\pi }}{2}\frac{C_{\text{dd}}(0)}{\Delta H},$$

$$k_{\mathrm{QI}}=\frac{\hslash \sqrt{\pi }}{2}\frac{C_{\text{dd}}(0)}{Q_r}\frac{C_{\text{ff}}(0)/C_{\text{dd}}(0)}{\Delta H}\text{,}$$

$$\Delta H=\hslash \sqrt{-\frac{{\ddot{C}}_{\text{dd}}(0)}{2C_{\text{dd}}(0)}}$$

$$\frac{\partial C_{\text{dd}}}{\partial {\eta }_k^{(\gamma )}}=0.$$

$$⟨{\mathbf{r}}^{(a)}|e^{-\beta \hat{H}}|{\mathbf{r}}^{(b)}⟩=⟨{\mathbf{r}}^{(a)}|(e^{-ϵ\hat{H}}{)}^P|{\mathbf{r}}^{(b)}⟩\text{.}$$

$$e^{-ϵ\hat{H}}=e^{-ϵ\hat{V}/2}{e}^{-ϵ\hat{T}}{e}^{-ϵ\hat{V}/2}+O\left({ϵ}^3\right)\text{.}$$

$$\mathrm{Tr}\left(e^{-ϵ\hat{H}}\right)=\mathrm{Tr}\left(e^{-ϵ{\hat{V}}_{\mathrm{TI}}/2}{e}^{-ϵ\hat{T}}{e}^{-ϵ{\hat{V}}_{\mathrm{TI}}/2}\right)+O\left({ϵ}^5\right),$$

$${\hat{V}}_{\mathrm{TI}}:=\hat{V}+\frac{1}{24}{ϵ}^2[\hat{V},[\hat{T},\hat{V}]]$$

$$e^{-ϵ\hat{H}}\approx e^{-ϵ{\hat{V}}_{\mathrm{TI}}/2}{e}^{-ϵ\hat{T}}{e}^{-ϵ{\hat{V}}_{\mathrm{TI}}/2}$$

$$e^{-ϵ\hat{H}}=e^{-ϵ{\hat{V}}_{\mathrm{e}}/6}{e}^{-ϵ\hat{T}/2}{e}^{-2ϵ{\hat{V}}_{\mathrm{m}}/3}{e}^{-ϵ\hat{T}/2}{e}^{-ϵ{\hat{V}}_{\mathrm{e}}/6}+O\left({ϵ}^5\right),$$

$${\hat{V}}_{\mathrm{e}}$$

$${\hat{V}}_{\mathrm{m}}$$

$$⟨{\mathbf{r}}^{(a)}|e^{-\beta \hat{H}}|{\mathbf{r}}^{(b)}⟩=\underset{P\to \infty }{\lim }C\int d{\mathbf{r}}^{(1)}\cdot \cdot \cdot d{\mathbf{r}}^{(P-1)}\exp \left[-\beta \tilde{\Phi }({\mathbf{r}}^{(a)},{\mathbf{r}}^{(1)},...,{\mathbf{r}}^{(P-1)},{\mathbf{r}}^{(b)})\right]\text{,}$$

$$\tilde{\Phi }$$

$$C$$

$$V_{\mathrm{eff}}^{(s)}:=V+{\left(\frac{\beta }{P}\right)}^2{d}_s{V}_{\mathrm{grad}},$$

$$V_{\mathrm{grad}}(\mathbf{r})={\hslash }^2||\nabla V(\mathbf{r})|{|}_{-}^2$$

$$d_s=\{\begin{array}{ll}0,& \text{LT splitting,}\\ 1/24,& \text{TI splitting,}\\ \alpha /6,& \text{SC splitting and }s\text{ even,}\\ (1-\alpha )/12,& \text{SC splitting and }s\text{ odd.}\end{array}$$

$$Q_r$$

$$Q_{r,P}$$

$${\rho }_r(\{{\mathbf{r}}^{(s)}\})$$

$$\Phi (\{{\mathbf{r}}^{(s)}\}):=\tilde{\Phi }({\mathbf{r}}^{(P)},{\mathbf{r}}^{(1)},...,{\mathbf{r}}^{(P)})=\frac{P}{2{\hslash }^2{\beta }^2}\sum _{s=1}^P||{\mathbf{r}}^{(s)}-{\mathbf{r}}^{(s-1)}|{|}_{+}^2+\frac{1}{P}\sum _{s=1}^P{w}_s{V}_{\mathrm{eff}}^{(s)}({\mathbf{r}}^{(s)})\text{.}$$

$$C_{\mathrm{dd},P}$$

$${\rho }^{‡}(\{{\mathbf{r}}^{(s)}\})$$

$${\overline{C}}_{\text{dd}}(\zeta ):=C_{\text{dd}}\left(-\frac{\mathrm{i}\zeta \hslash }{2}\right)$$

$${\ddot{C}}_{\text{dd}}(t)=-\frac{4}{{\hslash }^2}{\frac{{\partial }^2}{\partial {\zeta }^2}{\overline{C}}_{\text{dd}}\left(\zeta \right)|}_{\zeta =2\mathrm{i}t/\hslash }\text{.}$$

Full text

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 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.Sen and Kohen (2010)

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,Ramazani (2013) others use an approximation for the propagator by treating it semiclassicallyOlsson, Siegbahn, and Warshel (2004) or by treating only one or two degrees of freedom quantum mechanically.Banks and Clary (2007) Another promising method is the ring polymer molecular dynamicsAllen et al. (2013) (RPMD), which can partially capture both quantum effects and classical recrossing. Finally, there are various quantum generalizations of the transition state theory. Among these so-called quantum transition state theoriesHansen and Andersen (1994); Pollak (2012) belongs the quantum instanton (QI) approximation to the rate constant,Miller et al. (2003) i.e., the method whose efficiency we attempt to increase in the present paper. The QI approximation is motivated by the semiclassical instanton theory Miller (1974); Chapman, Garrett, and Miller (1975); Miller, Schwartz, and Tromp (1983); Ceotto (2012) 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.Zhao, Yamamoto, and Miller (2004); Ceotto and Miller (2004); Wang, Feng, and Zhao (2007); Wang and Zhao (2009, 2011, 2012)

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.Yamamoto and Miller (2004) As for KIE’s, the problem can be simplified further by using thermodynamic integration.Vaníček et al. (2005) 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.Jang, Jang, and Voth (2001); Yamamoto (2005); Pérez and Tuckerman (2011); Buchowiecki and Vaníček (2013); Marsalek et al. (2014); Engel et al. (2015) The second approach uses improved estimators with lower statistical errors, which permit shortening the MC simulation.Yang, Yamamoto, and Miller (2006); Vaníček and Miller (2007) 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 $\mathrm{\cdot H_{\alpha}+H_{\beta}H_{\gamma}\rightarrow H_{\alpha}H_{\beta}+\cdot H_{\gamma}}$ rearrangement, for which we also evaluate the resulting gain in computational efficiency, and the reaction $\mathrm{CH_{4}+\cdot H\rightleftharpoons\cdot CH_{3}+H_{2}}$, a process whose KIE’s were studied in detail both experimentally and theoretically, with classical TST, several of its corrected versions,Pu and Truhlar (2002); Schatz, Wagner, and Dunning (1984) reduced dimensionality quantum dynamics,Kerkeni and Clary (2004) and RPMD.Li et al. (2013)

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 the central 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 1.

II Quantum instanton formalism

The QI approximation for the thermal rate constant $k(T)$ can be derived from the exact Miller-Schwartz-Tromp formula,Miller, Schwartz, and Tromp (1983)

[TABLE]

expressing the product of the rate constant with the reaction partition function $Q_{r}$ as the time integral of the flux-flux correlation function

[TABLE]

where

[TABLE]

is the symmetrized correlation function of operators $\hat{A}$ and $\hat{B}$ at temperature $T=1/(k_{B}\beta)$ and time $t$,

[TABLE]

is the operator of flux through dividing surface (DS) $\ \gamma\in\{a,b\}$, $\mathbf{r}$ is the position vector in the $ND$-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\geq 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 $\xi_{\gamma}(\mathbf{r})=0$. In addition, $\xi_{\gamma}$ are chosen so that $\xi_{\gamma}(\mathbf{r})>0$ for $\mathbf{r}$ in the product region and $\xi_{\gamma}(\mathbf{r})<0$ in the reactant region.

The QI approximation can be derived by applying the steepest descent approximation to Eq. (1). Vaníček et al. (2005); Ceotto and Miller First, one multiplies and divides the integrand of Eq. (1) by the so-called delta-delta correlation function

[TABLE]

where $\Delta_{\gamma}$ is the normalized delta function

[TABLE]

and $||\cdot||_{-}$ is the norm of a covariant vector (see Table 1). Then one assumes that $C_{\text{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_{\text{ff}}(t)/C_{\text{dd}}(t)$ remains approximately constant and given by $C_{\text{ff}}(0)/C_{\text{dd}}(0)$. One can therefore evaluate the time integral in Eq. (1) with the steepest descent approximation,

[TABLE]

obtaining the QI expression for the rate constant

[TABLE]

where

[TABLE]

is a certain energy variance. For reasons that will become clear below, we keep $C_{\text{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_{\text{dd}}(0)$ be a saddle point with respect to $\xi_{a}$ and $\xi_{b}$;Miller et al. (2003) if $\xi_{\gamma}$ are controlled by a set of parameters $\{\eta_{k}^{(\gamma)}\},$ the stationarity condition becomes

[TABLE]

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,Feynman and Hibbs (1965); Chandler and Wolynes (1981) 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 Subsec. III.1, 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 Subsec. III.2 we obtain general PI expressions for $Q_{r}$ and $C_{\mathrm{dd}}(0)$. In Subsecs. III.3 and III.4, 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_{\mathrm{dd}}(0)$; in Subsec. III.5 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 the Appendix.

III.1 Lie-Trotter, Takahashi-Imada, and Suzuki factorizations of the

imaginary-time path integral propagator

The coordinate matrix element of the Boltzmann operator at temperature $T=1/(k_{B}\beta)$ can be rewritten as a matrix element of the product of $P\in\mathbb{N}$ Boltzmann operators at a higher temperature inversely proportional to the parameter $\epsilon:=\beta/P$:

[TABLE]

We next consider three possible high-temperature factorizations of the Boltzmann operator:

1. The symmetrized version of the Lie-Trotter factorization:

[TABLE]

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 Takahashi-Imada (TI) factorization:Takahashi and Imada (1984)

[TABLE]

where

[TABLE]

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 $\operatorname{Tr}$ in Eq. (13), and applies the resulting factorization

[TABLE]

to off-diagonal elements, which are required for PI representations of $C_{\text{dd}}(0)$ and $C_{\text{ff}}(0)$, one obtains only second-order 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. (13).

3. The fourth-order Suzuki-Chin (SC) factorization (Ref. Chin, 1997, motivated by Ref. Suzuki, 1995):

[TABLE]

where

[TABLE]

are the “endpoint” and “midpoint” effective one-particle potentials. The dimensionless parameter $\alpha$ can assume an arbitrary value, but evidence in the literatureJang, Jang, and Voth (2001); Pérez and Tuckerman (2011) suggests that $\alpha=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 $\epsilon=\beta/P$ with $\epsilon=2\beta/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

[TABLE]

where the effective potential $\tilde{\Phi}$ and prefactor $C$ are defined as

[TABLE]

In the expression for $\tilde{\Phi}$, we use the notation $\mathbf{r}^{(P)}:=\mathbf{r}^{(b)}$, $\mathbf{r}^{(0)}:=\mathbf{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$, $||\cdot||_{+}$ is the norm of a contravariant vector (see Table 1), and $V_{\mathrm{eff}}^{(s)}$ is the effective one-particle potential,

[TABLE]

where

[TABLE]

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:

[TABLE]

The weights $\tilde{w}_{s}$ in the sum over effective one-particle potentials also depend on the splitting: for endpoint $s$ (i.e., $s=0,P$) these weights are $\tilde{w}_{s}=1/2$ for the LT and TI splittings, and $\tilde{w}_{s}=1/3$ for the SC splitting; for other values of $s$, $\tilde{w}_{s}=1$ for the LT and TI splittings, whereas for the SC splitting, $\tilde{w}_{s}=4/3$ for odd $s$ and $\tilde{w}_{s}=2/3$ for even $s$. Expression (19) becomes exact as $P$ goes to infinity.

III.2 Path integral representation of the partition function and

delta-delta correlation function

From Eq. (19) it is straightforward to obtain the PI representation $Q_{r,P}$ of the reactant partition function $Q_{r}$; in particular,

[TABLE]

(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 $\int\mathit{d}\{\mathbf{r}^{(s)}\}$ we mean integration over all $\mathbf{r}^{(s)}$, $s\in\{1,...,P\}$; $\rho_{r}(\{\mathbf{r}^{(s)}\})$ can be regarded as the thermal distribution of $\{\mathbf{r}^{(s)}\}$ of the new system; $\Phi$ is the closed-loop version of $\tilde{\Phi}$, i.e.,

[TABLE]

From now on we will always consider closed loops such that $\mathbf{r}^{(0)}=\mathbf{r}^{(P)}$. The difference between the new weights $w_{s}$ and the old weights $\tilde{w}_{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,\text{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_{\text{dd}}(0)$ we have, analogously,

[TABLE]

(we shall omit the time argument of $C_{\text{dd}}$ and $C_{\text{ff}}$ if it equals 0). Note that $\rho^{\ddagger}(\{\mathbf{r}^{(s)}\})$ differs from $\rho_{r}(\{\mathbf{r}^{(s)}\})$ by the two delta constraints imposed on $\mathbf{r}^{(P/2)}$ and $\mathbf{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 $\rho_{r}$ and $\rho^{\ddagger}$. 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 $\rho_{r}$ and $\rho^{\ddagger}$ obtained with different Boltzmann operator splittings. While the use of the LT splitting only requires potential energy for each $\mathbf{r}^{(s)}$, the SC splitting with $0<\alpha<1$ and the TI factorization also require the gradient of energy for each $\mathbf{r}^{(s)}$, and the SC splittings with $\alpha=0$ and $\alpha=1$ require gradients for $\mathbf{r}^{(s)}$ with $s$ odd and even, respectively.

III.3 Estimators for constrained quantities

Within the PI formalism both the energy spread $\Delta H$ and the flux factor $C_{\mathrm{ff}}/C_{\mathrm{dd}}$ can be expressed as thermodynamic averages over the ensemble whose configurations are weighted by $\rho^{\ddagger}(\{\mathbf{r}^{(s)}\})$.Yamamoto and Miller (2004) In order to obtain the PI representation of $\ddot{C}_{\text{dd}}(0)$ and $\Delta H^{2}$, it is convenient to perform the Wick rotation and define a new function

[TABLE]

of a complex argument $\zeta$. Supposing that $C_{\text{dd}}(t)$ is analytic,

[TABLE]

The PI representation of $\overline{C}_{\text{dd}}(\zeta)$ is

[TABLE]

with

[TABLE]

prefactor

[TABLE]

and two “partial” effective potentials

[TABLE]

where $\tilde{\tilde{w}}_{s}=\tilde{w}_{s}$ for all $s$ except for $s=P/2$, for which $\tilde{\tilde{w}}_{P/2}=\tilde{w}_{P}=\tilde{w}_{0}$. The effective potentials $\tilde{\Phi}^{+}$ and $\tilde{\Phi}^{-}$ in Eq. (33) are obtained in a similar manner as $\tilde{\Phi}$ was obtained in Eq. (19). The difference is that instead of the matrix element of the Boltzmann operator $\exp(-\beta\hat{H})$ one considers an element of $\exp(-\beta^{+}\hat{H}/2)$ or $\exp(-\beta^{-}\hat{H}/2)$, and the exponential operators are discretized into $P/2$ rather than $P$ parts. As a result, expressions (37)-(38) for $\tilde{\Phi}^{+}$ and $\tilde{\Phi}^{-}$ can be obtained from the one for $\tilde{\Phi}$ [Eq. (20)] if $\beta$ is replaced with $\beta^{+}/2$ and $\beta^{-}/2$, respectively, and $P$ is replaced with $P/2$. After differentiating expression (33) with respect to $\zeta$, using Eq. (32) to go from $\overline{C}_{\text{dd},P}^{\prime\prime}(\zeta)$ back to $\ddot{C}_{\text{dd}}(t)$, and noting that $\mathit{d}\zeta=\mathit{d}\beta^{+}=-\mathit{d}\beta^{-}$, one obtains

[TABLE]

with

[TABLE]

After the substitution of expressions (29) and (39) for $C_{\mathrm{dd},P}$ and $\ddot{C}_{\mathrm{dd},P}$ into the definition (9) of $\Delta H^{2}$, the estimator for $\Delta H^{2}$ takes the form

[TABLE]

if $\rho^{\ddagger}(\{\mathbf{r}^{(s)}\})$ 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_{\mathrm{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,Yamamoto and Miller (2004)

[TABLE]

The ratio $C_{\text{ff}}/C_{\text{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 as

[TABLE]

Combining $\hat{F}_{\gamma}$ with the PI representation of the Boltzmann operator, one obtainsYamamoto and Miller (2004)

[TABLE]

where $f_{\mathrm{v}}$ is the so-called velocity factor,

[TABLE]

$\mathbf{r}_{a}=\mathbf{r}^{(P/2)}$, $\mathbf{r}_{b}=\mathbf{r}^{(P)}$, $\langle\cdot,\cdot\rangle_{-}$ is the inner product of two covariant vectors, and $\langle\cdot,\cdot,\cdot\rangle_{--}$ the matrix product of a covariant matrix with two covariant vectors (see Table 1). Taking the ratio of PI representations (46) and (29) of $C_{\text{ff}}$ and $C_{\text{dd}}$ immediately yields the estimator for the ratio $C_{\text{ff}}/C_{\text{dd}}$:

[TABLE]

The thermodynamic estimator takes the formYamamoto and Miller (2004)

[TABLE]

where $\langle\cdot,\cdot\rangle_{0}$ is the inner product of a covariant and contravariant vectors (see Table 1), and we defined

[TABLE]

[TABLE]

III.4 Thermodynamic integration with respect to mass

The last ingredient needed for evaluating the QI rate constant (8) is the ratio $C_{\text{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,Vaníček et al. (2005) which is easy to understand from the explicit QI expression for the KIE,

[TABLE]

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 $C_{\text{dd}}^{(A)}/C_{\text{dd}}^{(B)}$ and $Q_{r}^{(A)}/Q_{r}^{(B)}$ by considering a continuous transformationVaníček et al. (2005); Zimmermann and Vaníček (2009); Perez and von Lilienfeld (2011); Ceriotti and Markland (2013) from $A$ to $B$ using a dimensionless parameter $\lambda\in[0,1]$ controlling atomic masses of the intermediate systems as

[TABLE]

Ratios $C_{\text{dd}}^{(A)}/C_{\text{dd}}^{(B)}$ and $Q_{r}^{(A)}/Q_{r}^{(B)}$ are rewritten in terms of their logarithmic derivatives, which are normalized quantities and, therefore, can be calculated with the Metropolis algorithm:

[TABLE]

[TABLE]

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_{r}^{(A)}/Q_{r}^{(B)}$, 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 linearly interpolating $m^{-1/2}_{i}$ instead of $m_{i}$),Ceriotti and Markland (2013); Marsalek et al. (2014) or introducing higher-order derivatives of $Q_{r}$ with respect to $\lambda$;Marsalek et al. (2014) if the ratio is close to unity it is also possible to eliminate the integration error altogether by using direct estimators for $Q_{r}^{(A)}/Q_{r}^{(B)}$.Cheng and Ceriotti (2014)

In the case of $d\mathrm{ln}C_{\text{dd}}/\mathit{d}\lambda$ needed in the QI, one needs to keep track of the possible change of $\xi_{\gamma}$ during the course of the integration,

[TABLE]

In Ref. Vaníček et al., 2005 the authors proposed to choose $\{\eta_{k}^{(\gamma)}(\lambda)\}$ that satisfy Eq. (10) at each $\lambda$ integration step, making the second term in Eq. (56) exactly zero and leaving only $\partial\mathrm{ln}C_{\text{dd}}/\partial\lambda$ to be considered. Here we take an alternative and more numerically stable approach: By introducing new accurate estimators for $\partial\mathrm{ln}C_{\text{dd}}/\partial\eta_{k}^{(\gamma)},$ we can avoid having to find the optimal values of $\eta_{k}^{(\gamma)}(\lambda)$ for all $\lambda$. Instead, we only find optimal $\eta_{k}^{(\gamma)}(\lambda)$ at the boundary points $\lambda\in\{0,1\}$, obtain other, not necessarily optimal, $\eta_{k}^{(\gamma)}(\lambda)$ by linear interpolation, and evaluate both terms of Eq. (56) for each $\lambda$.

The estimators for $\partial\mathrm{ln}Q_{r}/\partial\lambda$ and $\partial\mathrm{ln}C_{\text{dd}}/\partial\lambda$ are

[TABLE]

where $F_{\mathrm{ds}}$ is the contribution that comes from differentiating the mass-dependent normalization factor in Eq. (6):

[TABLE]

Here $\nabla_{i}$ is the gradient with respect to coordinates of particle $i$ (see Table 1).

Direct evaluation of Eqs. (57) and (58) yields the thermodynamic estimatorsVaníček et al. (2005)

[TABLE]

Derivation of the estimator for $\partial\mathrm{ln}C_{\text{dd}}/\partial\eta_{k}^{(\gamma)}$ involves a rather tedious algebra and is therefore presented in Appendix B; the result is

[TABLE]

where $\nabla^{(\gamma)}$ is the gradient with respect to $\mathbf{r}_{\gamma}$ and

[TABLE]

is the term associated with the change of configuration space volume satisfying the constraint. Obtaining the thermodynamic estimator for $B^{k(\gamma)}$ is straightforward and yields

[TABLE]

III.5 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,Herman, Bruskin, and Berne (1982); Parrinello and Rahman (1984) which are motivated by the virial theorem of classical mechanics and which can be derivedPredescu and Doll (2002); Yamamoto (2005) most simply by applying a coordinate transformation before the differentiation.

Two of the five virial estimators used in this work, namely the estimators for $\partial\mathrm{ln}Q_{r}/\partial\lambda$ and $\Delta H^{2}$ had been proposed previously;Vaníček and Miller (2007); Yang, Yamamoto, and Miller (2006) the former, however, had not been used in combination with the SC factorization. To derive the centroid virial estimator for $\partial\mathrm{ln}Q_{r}/\partial\lambda$, let us choose an arbitrary bead $u$ and rewrite $Q_{r}$ in terms of the coordinates

[TABLE]

where $\{m_{i}^{\prime}\}$ are a set of parameters with dimensionality of mass. One then substitutes the new $C$ and $\Phi$ resulting from the transformation of coordinates into Eq. (57), and, finally, sets $\{m_{i}^{\prime}\}=\{m_{i}\}$ and transforms back to initial coordinates. This procedure yields an improved virial estimator,

[TABLE]

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\in\{1,...,P\}$ the centroid virial estimator is obtained,

[TABLE]

where

[TABLE]

is the centroid coordinate. 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. Vaníček and Miller, 2007, which was not considered in this work since its statistical error is larger then the error of its centroid counterpart.

For $\Delta H^{2}$, one startsYang, Yamamoto, and Miller (2006) by rewriting Eq. (33) using the coordinates

[TABLE]

where $\check{\mathbf{r}}^{(s)}$ is the reference point given by

[TABLE]

The kinetic parts of $\tilde{\Phi}^{\pm}$ are rewritten in the new coordinates; e.g., for $\tilde{\Phi}^{+}$, one uses the relation

[TABLE]

By substituting transformed $\tilde{\Phi}^{\pm}$ and $\overline{C}$ into Eqs. (40) and (41), one obtains the desired $G$ and $F$ terms of the virial estimator:

[TABLE]

[TABLE]

where $\langle\cdot,\cdot,\cdot\rangle_{00}$ is the matrix product of a covariant matrix with two contravariant vectors (see Table 1).

Now let us turn to the derivation of the new estimators promised in the Introduction. In particular, we propose new virial estimators for $\partial\mathrm{ln}C_{\mathrm{dd}}/\partial\lambda$, $C_{\mathrm{ff}}/C_{\mathrm{dd}}$, and $\partial\mathrm{ln}C_{\text{dd}}/\partial\eta_{k}^{(\gamma)}$. For $\partial\mathrm{ln}C_{\text{dd}}/\partial\lambda$ we use a coordinate rescaling

[TABLE]

which is similar to Eq. (69) and yields the virial estimator

[TABLE]

For $C_{\text{ff}}/C_{\text{dd}}$, we introduce new coordinates

[TABLE]

and employ the identity

[TABLE]

Rewriting $\tilde{\Phi}^{\pm}$ in terms of $\{\mathbf{x}^{(1)},$ $...,$ $\mathbf{x}^{(P/2-1)},$ $\mathbf{r}^{(P/2)},$ $\mathbf{x}^{(P/2+1)},$ $...,$ $\mathbf{x}^{(P-1)},$ $\mathbf{r}^{(P)}\}$ and inserting them into Eq. (47) leads to the virial estimator

[TABLE]

where we introduced coefficients

[TABLE]

Using the same rescaling as for $f_{\mathrm{v,v}}$, we can also derive the virial estimator for $\partial\mathrm{ln}C_{\text{dd}}/\partial\eta_{k}^{(\gamma)}$,

[TABLE]

where $\mathbf{r}_{\gamma}^{(P/2)}$ stands for $\mathbf{\mathbf{r}}^{(P)}$ if $\gamma=a$ and for $\mathbf{\mathbf{r}}^{(P/2)}$ if $\gamma=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_{\mathrm{th}}^{k(\gamma)}$ 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, $\langle\mathbf{w},\nabla V\rangle_{0}$ and $\langle\mathbf{w},\partial^{2}V/\partial\mathbf{r}^{2},\mathbf{w}\rangle_{00}$ are first and second derivatives of $V$ in the direction of $\mathbf{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,\mathrm{cv}}$, one per unconstrained bead for $F_{\mathrm{cv}}^{\ddagger}$, two per unconstrained bead for $(G_{\mathrm{v}}+F_{\mathrm{v}}^{2})/2$, and three for $f_{\mathrm{v,v}}$. Calculating $B_{\mathrm{th}}^{k(\gamma)}$ will require exactly one potential evaluation and calculating $B_{\mathrm{v}}^{k(\gamma)}$ will require $P-1$ evaluations unless it is computed at the same time as $f_{\mathrm{v,v}}$ (in this case it would require just one extra potential evaluation, other numerical ingredients being shared with $f_{\mathrm{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. Vaníček and Miller, 2007 and Yamamoto, 2005 used finite differences with respect to mass and $\beta$, 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:

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 $\xi_{\gamma}(\mathbf{r})$ of the form

[TABLE]

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.,

[TABLE]

where $R_{XY}$ is the distance between X and Y. Optimal values of $\eta_{\gamma}^{\ddagger}$ were found by running test simulations to find the sign of $\partial\mathrm{ln}C_{\text{dd}}/\partial\eta_{\gamma}^{\ddagger}$ at different values of $(\eta_{a}^{\ddagger},\eta_{b}^{\ddagger})$.

3. Run simulations at different values of $\lambda$ in order to obtain the corresponding logarithmic derivatives of $Q_{r}$ and $C_{\text{dd}}$, as well as $C_{\text{ff}}(0)/C_{\text{dd}}(0)$ and $\Delta H$ for $\lambda=0$ and $\lambda=1$, then evaluate Eqs. (54) and (55) using, e.g., Simpson’s rule. For many systems $d\mathrm{ln}C_{\text{dd}}/\mathit{d}\lambda$ and $\partial\mathrm{ln}Q_{r}/\partial\lambda$ are quite smooth functions and nine intermediate points were sufficient to accurately evaluate the thermodynamic integrals (i.e., the discretization error of the $\lambda$ integral was smaller than the already small statistical error). After this, evaluating the KIE using Eq. (52) is straightforward.

For each value of $\lambda$ one has to run two MC simulations in $\{\mathbf{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 $\delta_{\text{sm}}$,

[TABLE]

In contrast with the approximation used in Ref. Yamamoto and Miller, 2004, the width $\sigma$ of our Gaussian $\delta_{m}$ does not depend on temperature or mass. The approximate constraint converges to the exact delta function as $\sigma\rightarrow 0$. Presence of $\delta_{\text{sm}}[\xi_{\gamma}(\mathbf{r})]$ can be easily simulated by adding an extra constraining potential to two of the slices. For MC sampling, we employed the staging algorithmSprik, Klein, and Chandler (1985a, b); Tuckerman et al. (1993) 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.

IV.1 $\mathbf{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 $\lambda$ steps.) To verify the improvements outlined in Sec. III we studied their influence on the behavior of the two main types of errors when applied to the model $\mathrm{\cdot H_{\alpha}+H_{\beta}H_{\gamma}\rightarrow H_{\alpha}H_{\beta}+\cdot H_{\gamma}}$ rearrangement using the BKMP2 potential energy surfaceBoothroyd et al. (1996) at the temperature of $200\,\text{K}$. The behavior of the logarithmic derivatives was studied on the KIE $\mathrm{\cdot H+H_{2}/\cdot D+D_{2}}$.

IV.1.1 Computational details

Statistically converged simulations (paralleled over 64 trajectories, $4\times 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,Flyvbjerg and Petersen (1989) 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 $\eta_{a}^{\ddagger}=\eta_{b}^{\ddagger}=0$ was quite satisfactory even at $T=200$ K (in this case $C_{\text{dd}}$ is stationary from symmetry considerations) for analyzing numerical behavior of $\partial\mathrm{ln}C_{\text{dd}}/\partial\lambda$, $\Delta H^{2}$ and $C_{\mathrm{ff}}/C_{\mathrm{dd}}$. For $\partial\mathrm{ln}C_{\text{dd}}/\partial\eta_{a}^{\ddagger}$, however, we used $\eta_{a}^{\ddagger}=-0.5$ and $\eta_{b}^{\ddagger}=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.

IV.1.2 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 $\partial\mathrm{ln}Q_{r}/\partial\lambda$ 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 $\Delta H$, $C_{\mathrm{ff}}/C_{\mathrm{dd}}$, $Q_{r}^{(B)}/Q_{r}^{(A)}$, and $C_{\text{dd}}^{(B)}/C_{\text{dd}}^{(A)}$ to $1\%$ discretization and statistical errors. 111Note that the $1\%$ error for $\Delta H$ translates into a $2\%$ error for $\Delta H^{2}$ and that $1\%$ relative error for $Q_{r}$ and $C_{\text{dd}}$ ratios translate into $0.01$ absolute error for $\partial\ln Q_{r}/\partial\lambda$ and $\partial\ln C_{\mathrm{dd}}/\partial\lambda$. As for $\partial\ln C_{\mathrm{dd}}/\partial\eta_{a}^{\ddagger}$, as will be shown later, when we calculate the KIE $\mathrm{\cdot H+H_{2}/\cdot D+D_{2}}$ at $T=200\,\text{K}$ with a properly optimized DS, $\partial\ln C_{\mathrm{dd}}/\partial\eta_{a}^{\ddagger}$ is integrated over an interval of the length $0.59$ a. u., implying that the target error should be $0.01/0.59(\text{a.\thinspace u.})^{-1}$. 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 $\lambda=1$ in addition to those for $\lambda=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_{\mathrm{dd}}$ had been performed using a single step trapezoidal rule (i.e., based just on the two boundary points $\lambda=0$ and $\lambda=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 product $\tilde{P}\sigma_{\tilde{P}}^{2}$, where $\tilde{P}$ is the value of the Trotter number that yields the target discretization error and $\sigma_{\tilde{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 Tab. IV.1.2, 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_{\mathrm{ff}}/C_{\mathrm{dd}}$ appears to be roughly independent of the splitting used; this is probably because the discretization errors of $C_{\mathrm{dd}}$ and $C_{\mathrm{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 $\partial\ln Q_{r}/\partial\lambda$ 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 KIE computed with the LT splitting is still larger than, for instance, those used in Ref. Yamamoto and Miller, 2004, where the authors obtained the final result by extrapolating to the $P\rightarrow\infty$ limit.222Since thermodynamic estimators were used in Ref. Yamamoto and Miller, 2004, reducing the discretization error directly using very large $P$ was not feasible—increasing $P$ not only decreased discretization error, but also increased the statistical error. Introducing virial estimators for each relevant quantity allows avoiding this issue because it permits improving convergence with respect to $P$ without encountering problems with statistical error.

It is also interesting to relate our results to those of Ref. Engel et al., 2015, where the authors compared efficiencies of the LT, TI, and fourth-order ChinChin and Chen (2002); Chin (2004) 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 Subsec. III.1; both factorizations are outperformed by the fourth-order Chin factorization, which is in agreement with the SC outperforming LT splitting in Tab. IV.1.2. 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 $\alpha=0$ or $\alpha=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 Tab. IV.1.2, 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;Jang, Jang, and Voth (2001); Yamamoto (2005); Pérez and Tuckerman (2011) this approach, however, is known to increase the statistical errors of the final result in high-dimensional systems.Ceriotti et al. (2012)