---
abstract: 'A general and rigorous methodology to compute the quantum equilibrium isotope effect is described. Unlike standard approaches, ours does not assume separability of rotational and vibrational motions and does not make the harmonic approximation for vibrations or rigid rotor approximation for the rotations. In particular, zero point energy and anharmonicity effects are described correctly quantum mechanically. The approach is based on the thermodynamic integration with respect to the mass of isotopes and on the Feynman path integral representation of the partition function. An efficient estimator for the derivative of free energy is used whose statistical error is independent of the number of imaginary time slices in the path integral, speeding up calculations by a factor of $\sim60$ at $500$ $\mathrm{K}$ and more at room temperature. We describe the implementation of the methodology in the molecular dynamics package Amber 10. The method is tested on three [\[]{}1,5\] sigmatropic hydrogen shift reactions. Because of the computational expense, we use ab initio potentials to evaluate the equilibrium isotope effects within the harmonic approximation, and then the path integral method together with semiempirical potentials to evaluate the anharmonicity corrections. Our calculations show that the anharmonicity effects amount up to $30\%$ of the symmetry reduced reaction free energy. The numerical results are compared with recent experiments of Doering and coworkers, confirming the accuracy of the most recent measurement on 2,4,6,7,9-pentamethyl-5-(5,5-$^{2}$H$_{2}$)methylene-11,11a-dihydro-12*H*-naphthacene as well as concerns about compromised accuracy, due to side reactions, of another measurement on 2-methyl-10-(10,10-$^{2}$H$_{2}$)methylenebicyclo[\[]{}4.4.0\]dec-1-ene.'
author:
- Tomáš Zimmermann
- Jiří Vaníček
title: Path integral evaluation of equilibrium isotope effects
---

Introduction\[sec:Introduction\]
================================

The equilibrium (thermodynamic) isotope effect (EIE) is defined as the effect of isotopic substitution on the equilibrium constant. Denoting an isotopolog with a lighter (heavier) isotope by a subscript *l* (*h*), the EIE is defined as the ratio of equilibrium constants$$\mathrm{EIE}=\frac{K_{l}}{K_{h}}=\frac{Q_{l}^{\left(p\right)}/Q_{l}^{\left(r\right)}}{Q_{h}^{\left(p\right)}/Q_{h}^{\left(r\right)}},\label{eq:EIE_def}$$ where $Q^{\left(r\right)}$ **and $Q^{\left(p\right)}$ are molecular partition functions per unit volume of reactant and product. We study a specific case of EIE - the equilibrium ratio of two isotopomers. In this case, the EIE is equal to the equilibrium constant of the isotopomerization reaction, $$\mathrm{EIE}=K_{\mathrm{eq}}=\frac{Q^{\left(p\right)}}{Q^{\left(r\right)}},\label{eq:EIE}$$ where superscripts *r* and *p* refer to reactant and product isotopomers, respectively.

Usually, equilibrium isotope effects are computed only approximately:[@Alston1996; @Anbar2005825; @Gawlita1994; @Hill20081939; @Hrovat1997; @Janak2003a; @Kolmodin2002; @Munoz-Caro2003; @Ruszczycky2009107; @Saunders2007; @Slaughter2000; @Smirnov2009; @Zeller2002] In particular, effects due to indistinguishability of particles and rotational and vibrational contributions to the EIE are treated separately. Furthermore, the vibrational motion is approximated by a simple harmonic oscillator and the rotational motion is approximated by a rigid rotor. In general, none of the contributions, not even the indistinguishability effects can be separated from the others.[@Pollock1987; @Boninsegni1995] However, at room temperature or above the nuclei can be accurately treated as distinguishable, and the indistinguishability effects can be almost exactly described by symmetry factors. On the other hand, the effective coupling between rotations and vibrations, anharmonicity of vibrations, and non-rigidity of rotations can in fact become more important at higher temperatures. For simplicity, from now on we denote these three effects together as “anharmonicity effects and the approximation that neglects them the “harmonic approximation (HA). In some cases the effects of anharmonicity of the Born-Oppenheimer potential surface on the value of EIE can be substantial.[@TorresL._jp001327p] Ishimoto *et al.* have shown that the isotope effect on certain barrier heights[^1] can even have opposite signs, when calculated taking anharmonicity effects into account and in the HA.[@Ishimoto2008]

Our goal is to describe rigorously equilibria at room temperature or above. Therefore, two approximations that we make are the Born-Oppenheimer approximation and the distinguishable particle approximation (we treat indistinguishability by appropriate symmetry factors). The error due to the Born-Oppenheimer approximation was studied for H/D EIE by Bardo and Wolfsberg,[@Bardo1976] and by Kleinman and Wolfsberg,[@Kleinman1973] and was shown to be of the order of 1 % in most studied cases. Since we assume that nuclei are point charges, the Born-Oppenheimer approximation implies that the potential energy surface is the same for the two isotopomers. The differences of Born-Oppenheimer surfaces due to differences of nuclear volume and quadrupole of isotopes can be important for heavy elements[@Bigeleisen1996; @Knyazev19992579; @Cowan1981], but these are not studied in this work.

Symmetry factors themselves result in an EIE equal to a rational ratio, which can be computed analytically. In order to separate the symmetry contributions from the mass contributions to the EIE, it is useful to introduce the reduced reaction free energy,$$\Delta F^{\mathrm{red}}=-k_{B}T\ln K_{\mathrm{eq}}^{\mathrm{red}}=-k_{B}T\ln\left(\frac{s^{(r)}Q^{(p)}}{s^{(p)}Q^{(r)}}\right),\label{eq:F_red}$$ where $s^{(r)}$ and $s^{(p)}$ are the symmetry factors discussed in more detail in Sec. \[sec:Examples\]. To include the effects of quantization of nuclear degrees of freedom beyond the HA rigorously we use the Feynman path integral representation (PI) of the partition function. The quantum reduced reaction free energy can then be computed by thermodynamic integration with respect to the mass of the isotopes. To compute the derivative of the free energy efficiently, we use a generalized virial estimator. The advantage of this estimator is that its statistical error does not increase with the number of imaginary time slices in the discretized path integral. As a consequence, converged results can be obtained in a significantly shorter simulation than with other estimators.

The ultimate goal would be to combine the path integral methodology with *ab initio* potentials. However, since millions of samples are required, the computational expense results in the following “compromise: First, the equilibrium isotope effects are computed using *ab initio* potentials, but as usual, within the HA. Then all anharmonicity corrections are computed using the PI methodology, but with semiempirical potentials. In other words, we take advantage of the higher accuracy of the *ab initio* potentials to compute the harmonic contribution to the EIE and then make an assumption that the anharmonicity effects are similar for *ab initio* and semiempirical potentials.

After describing theoretical features of the method, we apply it to hydrocarbons used in experimental measurements of isotope effects on [\[]{}1,5\] sigmatropic hydrogen transfer reactions. Two of them were recently used by Doering *et al.* [@Doering2006; @Doering2007] who reported equilibrium ratios of their isotopomers. This allows us to validate our calculations as well as to discuss the apparent discrepancy in measurements of Doering *et al.* from a theoretical point of view.

The outline of the paper is as follows: In Sec. \[sec:Method\], we describe a rigorous quantum-mechanical methodology to compute the EIE. Section \[sec:Examples\] presents the [\[]{}1,5\] sigmatropic hydrogen shift reactions on which we test the methodology, explains how *ab initio* methods can be combined with the PI to compute the EIE, and discusses in detail symmetry effects in these reactions. Section \[sec:Computational-details\] explains the implementation of the method in Amber 10 and describes details of calculations and error analysis of our PIMD simulations. Results of calculations are presented and compared with experiments in Sec. \[sec:Results\]. Section \[sec:Discussion-and-conclusions\] concludes the paper.

The methodology\[sec:Method\]
=============================

Thermodynamic integration {#thermodynamic-integration .unnumbered}
-------------------------

EIE can be calculated by a procedure of thermodynamic integration[@Chandler:1987] with respect to the mass. This method takes advantage of the relationship$$\mathrm{EIE}=\frac{Q^{(p)}}{Q^{(r)}}=\exp\left[-\beta\int_{0}^{1}d\lambda\frac{dF\left(\lambda\right)}{d\lambda}\right],\label{eq:TI}$$ where $F=-\log Q/\beta$ is the (quantum) free energy and $\lambda$ is a parameter which provides a smooth transition between isotopomers $r$ and $p$. This can be accomplished, e.g., by linear interpolation of masses of all atoms in a molecule according to the equation $$m_{i}\left(\lambda\right)=\left(1-\lambda\right)m_{i}^{(r)}+\lambda m_{i}^{(p)}.\label{eq:lambda_mass}$$ In contrast to the partition function itself, the integrand of Eq. $$\frac{dF\left(\lambda\right)}{d\lambda}=-\frac{1}{\beta}\frac{d\log Q(\lambda)}{d\lambda}=-\frac{1}{\beta}\frac{dQ(\lambda)/d\lambda}{Q(\lambda)}\label{eq:log_derivative}$$ is a thermodynamic average and therefore can be computed by either Monte Carlo or molecular dynamics simulations.

Path integral approach {#path-integral-approach .unnumbered}
----------------------

Classically, the EIE is trivial and Eq. can be evaluated analytically. When quantum effects are important, this simplification is not possible. To describe quantum thermodynamic effects rigorously, one can use the path integral formulation of quantum mechanics.[@Feynman1965] In the path integral formalism, thermodynamic properties are calculated exploiting the correspondence between matrix elements of the Boltzmann operator and the quantum propagator in imaginary time.[@Feynman1965; @Pollock1987] In the last decades, path integrals proved to be very useful in many areas of quantum chemistry, most recently in calculations of heat capacities,[@Predescu2003] rate constants,[@yamamoto:044106] kinetic isotope effects,[@vanicek:114309] or diffusion coefficients.[@wang:114708]

Let $N$  be the number of atoms, $D$ the number of spatial dimensions, and $P$ the number of imaginary time slices in the discretized PI ($P=1$ gives classical mechanics, $P\rightarrow\infty$ gives quantum mechanics). Then the PI representation of the partition function $Q$ in the Born-Oppenheimer approximation is $$\begin{aligned}
Q & \simeq V^{-1}C\int d\mathbf{r}^{\left(0\right)}\cdots\int d\mathbf{r}^{\left(P-1\right)}\exp\left[-\beta\Phi\left(\left\{ \mathbf{r}^{\left(s\right)}\right\} \right)\right],\label{eq:pi_qr}\\
C & \equiv\left(\frac{P}{2\pi\hbar^{2}\beta}\right)^{NPD/2}\prod_{i=1}^{N}m_{i}^{PD/2}.\label{eq:pi_prefactor}\end{aligned}$$ where $\mathbf{r}^{\left(s\right)}\equiv\left(\mathbf{r}_{1}^{\left(s\right)},\mathbf{r}_{2}^{\left(s\right)},\ldots,\mathbf{r}_{N}^{\left(s\right)}\right)$ is the set of Cartesian coordinates associated with the $s$$^{\text{th}}$ time slice, and $\Phi\left(\left\{ \mathbf{r}^{\left(s\right)}\right\} \right)$ is the effective potential $$\begin{split}
\Phi\left(\left\{ \mathbf{r}^{\left(s\right)}\right\} \right)=\frac{P}{2\hbar^{2}\beta^{2}}\sum_{i=1}^{N}m_{i}\sum_{s=0}^{P-1}\left(\mathbf{r}_{i}^{\left(s\right)}-\mathbf{r}_{i}^{\left(s+1\right)}\right)^{2}\\+\frac{1}{P}\sum_{s=0}^{P-1}V\left(\mathbf{r}^{\left(s\right)}\right).\end{split}\label{eq:pi_effpotential}$$ The $P$ particles representing each nucleus in $P$ different imaginary time slices are called “beads.” Each bead interacts via harmonic potential with the two beads representing the same nucleus in adjacent time slices and via potential $V\left(\mathbf{r}^{\left(s\right)}\right)$ attenuated by factor $1/P$ with beads representing other nuclei in the same imaginary time slice. By straightforward differentiation of Eq. we obtain the so-called thermodynamic estimator (TE),[@vanicek:054108] $$\begin{split}
\frac{dF\left(\lambda\right)}{d\lambda}\simeq-\frac{1}{\beta}\sum_{i=1}^{N}\frac{dm_{i}}{d\lambda}\Bigg\langle \frac{DP}{2m_{i}}\\-\frac{P}{2\hbar^{2}\beta}\sum_{s=0}^{P-1}\left(\mathbf{r}_{i}^{\left(s\right)}-\mathbf{r}_{i}^{\left(s+1\right)}\right)^{2}\Bigg\rangle.\end{split}\label{eq:TE_EST}$$ A problem with expression is that its statistical error grows with the number of time slices. A similar behavior is a well known property of the thermodynamic estimator for energy,[@herman:5150] where the problem is caused by the kinetic part of energy.

![Convergence of the generalized virial estimator (GVE) and the thermodynamic estimator (TE) as a function of the number $P$ of imaginary time slices in the path integral. All results obtained by 1 $\mathrm{ns}$ long simulations (with the time step 0.05 $\mathrm{fs}$) of compound **1-5,5,5-*d*$_{\text{3}}$** ($\lambda=0$) using GAFF force field, normal mode PIMD, and Nosé-Hoover chains of thermostats. \[fig:GVE\_TE\_values\]](Fig1_jcp_300dpi)

Generalized virial estimator  {#generalized-virial-estimator .unnumbered}
-----------------------------

![Root mean square errors (RMSEs) of the generalized virial estimator (GVE) and the thermodynamic estimator (TE) as a function of the number $P$ of imaginary time slices in the path integral. All results obtained by 1 $\mathrm{ns}$ long simulations (with the time step 0.05 $\mathrm{fs}$) of compound **1-5,5,5-*d*$_{\text{3}}$** ($\lambda=0$) using GAFF force field, normal mode PIMD, and Nosé-Hoover chains of thermostats. Note that the RMSE of the GVE is not only non-increasing (as expected from theory), but in fact decreases slightly with increasing $P$, which is due to the decrease of correlation length.\[fig:GVE\_TE\_errors\]](Fig2_jcp_300dpi)

The growth of statistical error of the thermodynamic estimator for energy is removed by expressing the estimator only in terms of the potential and its derivatives using the virial theorem.[@herman:5150] In our case, a similar improvement can be accomplished, if a coordinate transformation, $$\mathbf{x}_{i}^{\left(s\right)}=m_{i}^{1/2}(\mathbf{r}_{i}^{\left(s\right)}-\mathbf{r}_{i}^{(C)}),\label{eq:mass_centroid_transf}$$ is introduced into Eq. prior to performing the derivative. Here, the “centroid” coordinate is defined as

$$\mathbf{r}_{i}^{(C)}=\frac{1}{P}\sum_{s=0}^{P-1}\mathbf{r}_{i}^{(s)}.\label{eq:centroid_coord}$$

In other words, first the “centroid” coordinate $\mathbf{r}_{i}^{\left(C\right)}$ is subtracted, and then the coordinates are mass-scaled. Resulting generalized virial estimator (GVE)[@vanicek:114309; @Vanicek_Dubna_2005] takes the form

$$\begin{split}
\frac{dF\left(\lambda\right)}{d\lambda}\simeq-\frac{1}{\beta}\sum_{i=1}^{N}\frac{dm_{i}/d\lambda}{m_{i}}\Bigg[\frac{D}{2}\\+\frac{\beta}{2P}\left\langle \sum_{s=0}^{P-1}\left(\mathbf{r}_{i}^{(s)}-\mathbf{r}_{i}^{(C)}\right)\cdot\frac{\partial V\left(\mathbf{r}^{(s)}\right)}{\partial\mathbf{r}_{i}^{(s)}}\right\rangle \Bigg].
\end{split}
\label{eq:gve}$$

Its primary advantage is that the root mean square error (RMSE) of the average, $\sigma_{\mathrm{\text{av}}}$, is approximately independent on the number of imaginary time slices $P$,$$\sigma_{\mathrm{\text{av}}}\approx O\left(P^{0}\tau_{\mathrm{c}}^{1/2}\tau_{\mathrm{sim}}^{-1/2}\right).\label{eq:sigma_virial}$$ In this equation $\tau_{\mathrm{sim}}$ denotes the length of the simulation and $\tau_{\mathrm{c}}$ is the correlation length.

The convergence of values and statistical errors of both estimators as a function of number $P$ of imaginary time slices for systems studied in this paper is discussed in Section \[sec:Computational-details\]. As expected, up to the statistical error they give the same values as can be seen in Fig. \[fig:GVE\_TE\_values\]. Nevertheless, when quantum effects are important, and a high value of $P$ must be used, the GVE is the preferred estimator since it has much smaller statistical error and therefore converges much faster than the TE (see Fig. \[fig:GVE\_TE\_errors\]).

Path integral molecular dynamics {#path-integral-molecular-dynamics .unnumbered}
--------------------------------

Thermodynamic average in Eq. can be evaluated efficiently using the path integral Monte Carlo (PIMC) or path integral molecular dynamics (PIMD). In PIMC, gradients of $V$ in Eq. result in additional calculations since the usual Metropolis Monte Carlo procedure for the random walk only requires the values of $V$. This additional cost can be, however, reduced either by less frequent sampling, or by using a trick in which the total derivative with respect to $\lambda$ (not the gradients!) is computed by finite difference.[@vanicek:114309; @Vanicek_Dubna_2005; @Predescu2003; @Predescu2004] In case of PIMD, the presence of gradients of $V$ in Eq. does not slow down the calculation since forces are already computed by a propagation algorithm. Although in principle, a PIMC algorithm for a specific problem can always be at least as efficient as a PIMD algorithm, in practice it is much easier to write a general PIMD algorithm and so PIMD is usually the algorithm used in general software packages. Since PIMD was implemented in Amber 9,[@amber9] we have implemented the methodology described above for computing EIEs into Amber 10.[@amber10] This implementation is what was used in calculations in the following sections. In PIMD, equation is augmented by fictitious classical momenta $\boldsymbol{\mathrm{p}}^{\left(s\right)}$,

$$Q\simeq V^{-1}\tilde{C}\int d\mathbf{p}^{\left(0\right)}\cdots\int d\mathbf{p}^{\left(P-1\right)}\int d\mathbf{r}^{\left(0\right)}\cdots\int d\mathbf{r}^{\left(P-1\right)}\exp\left[-\beta\left({\displaystyle \sum_{s=0}^{P-1}}\frac{\left(\mathbf{p}^{\left(s\right)}\right)^{2}}{2m_{s}}+\Phi\left(\left\{ \mathbf{r}^{\left(s\right)}\right\} \right)\right)\right].\label{eq:pi_qp}$$

The normalization prefactor $\tilde{C}$ is chosen in such a way that the original prefactor $C$ in Eq. is reproduced when the momentum integrals in Eq. are evaluated analytically. The partition function is formally equivalent to the partition function of a classical system of cyclic polyatomic molecules with harmonic bonds. Each such molecule represents an individual atom in the original molecule and interacts with molecules representing other atoms via a potential derived from the potential of the original molecule.[@Chandler1981] This system can be studied directly using well developed methods of the classical molecular dynamics.

Low- and high-temperature limits {#low--and-high-temperature-limits .unnumbered}
--------------------------------

It is useful to see how the general path integral expressions , , and behave in the low and high temperature limits. As temperature decreases, the difference of the reduced free energies of isotopomers approaches the difference of their zero point energies (ZPEs). Therefore, still assuming that indistinguishability is correctly described by symmetry factors, the low temperature limit of the EIE is equal to $$\mathrm{EIE_{low\mathit{T}}}=\frac{s^{\left(r\right)}}{s^{\left(p\right)}}\exp\left[-\beta\left(\varepsilon^{\left(r\right)}-\varepsilon^{\left(p\right)}\right)\right],\label{eq:EIE_low_temp}$$ where $\varepsilon$ denotes the ZPE.

At high temperature, the system approaches its classical limit. In this limit, we can set the number of imaginary time slices $P=1$, and the EIE can be computed analytically using partition function , which becomes$$\begin{split}
Q_{\mathrm{class}}\simeq {\displaystyle V^{-1}\left(\frac{1}{2\pi\hbar^{2}\beta}\right)^{ND/2}\prod_{i=1}^{N}m_{i}^{D/2}}\\\times\int d^{N}\mathbf{r}\exp\left[-\beta V\left(\mathbf{r}_{1},\ldots,\mathbf{r}_{N}\right)\right].\end{split}\label{eq:q_class}$$ Again, for isotopomers the mass dependent factors cancel out upon substitution into Eq. , along with all other terms except for the integrals. Although the Born-Oppenheimer potential surface is the same for all isotopomers, the integrals do not cancel out since the integration is not performed over the whole configuration space, but only over parts attributed to reactants or products. Particles are treated as distinguishable in the classical limit and the volumes of the configuration space corresponding to reactants and products are generally different. After the reduction, EIE becomes exactly equal to the ratio of the symmetry factors,

$$\mathrm{EIE_{\mathrm{high}\mathit{T}}}=\frac{s^{\left(r\right)}}{s^{\left(p\right)}}.\label{eq:EIE_class}$$

Equilibrium isotope effects in three [\[]{}1,5\] sigmatropic hydrogen shift reactions\[sec:Examples\]
=====================================================================================================

We examine the EIE for four related compounds. The parent compound (3*Z*)-penta-1,3-diene (compound **1**) is the simplest molecule to model the [\[]{}1,5\] sigmatropic hydrogen shift reaction. Two of its isotopologs, tri-deuterated (3*Z*)-(5,5,5-$^{\text{2}}$H$_{\text{3}}$)penta-1,3-diene (**1-5,5,5-*d*$_{\text{3}}$**) and di-deuterated (3*Z*)-(1,1-$^{\text{2}}$H$_{2}$)penta-1,3-diene (**1-1,1-*d*$_{\text{2}}$**) (see Fig. \[fig:equlibrium\_pent\]) were used by Roth and König to measure an unusually high value of the kinetic isotope effect (KIE) on the [\[]{}1,5\] hydrogen shift reaction with respect to the substitution of hydrogen by deuterium.[@Roth1966] This result pointed to a significant role of tunneling in [\[]{}1,5\] sigmatropic hydrogen transfer reactions. Subsequently, much theoretical research was devoted to the study of this reaction. Here we calculate final equilibrium ratios of products of this reaction, which, to our knowledge, were not theoretically predicted so far.

![Equilibrium of the [\[]{}1,5\] hydrogen shift reaction in (3*Z*)-(5,5,5-$^{\text{2}}$H$_{\text{3}}$)penta-1,3-diene (**1-5,5,5-*d*$_{\text{3}}$**) and in (3*Z*)-(1,1-$^{\text{2}}$H$_{2}$)penta-1,3-diene (**1-1,1-*d*$_{\text{2}}$**). If all contributions except for those due to symmetry factors $s_{a}$ and $s_{b}$ were neglected, one would obtain approximate equilibrium constants $K_{1}=3$ and $K_{2}=2$ in both cases.\[fig:equlibrium\_pent\]](Fig3_jcp_300dpi)

Two other compounds, 2-methyl-10-(10,10-$^{\text{2}}$H$_{2}$)methylenebicyclo[\[]{}4.4.0\]dec-1-ene (**2-1,1-*d*$_{\text{2}}$**) and 2,4,6,7,9-pentamethyl-5-(5,5-$^{\text{2}}$H$_{2}$)methylene-11,11a-dihydro-12*H*-naphthacene (**3-1,1-*d*$_{\text{2}}$**) (see Fig. \[fig:equilibria\_cyclic\]), were recently used by Doering *et al.* to confirm and possibly refine the experimental value of the KIE on the [\[]{}1,5\] hydrogen shift.[@Doering2006; @Doering2007] In contrast to (3*Z*)-penta-1,3-diene **(**compound **1**), where the s-*trans* conformer incompetent of the [\[]{}1,5\] hydrogen shift is the most stable, pentadiene moiety in compounds 2 and 3 is locked in the s-*cis* conformation. This not only increases the reaction rate, but also rules out the (very small) effect of the EIE on the KIE due to the shift in s-*cis*/s-*trans* equilibrium. For both molecules Doering *et al.* reported final equilibrium ratios of isotopomers. Despite the similarity of compounds **2** and **3**, these ratios are qualitatively different. Indeed, one motivation for measuring EIE in compound **3** was that Doering *et al.* suspected that in the case of **2-1,1-*d*$_{\text{2}}$** the equilibrium ratio might be modified by unwanted side reactions, mainly dimerizations. One of our goals is to elucidate this discrepancy from the theoretical point of view.

![image](Fig4_jcp_300dpi)

As can be seen in Figs. \[fig:equlibrium\_pent\] and \[fig:equilibria\_cyclic\], the final equilibrium of the [\[]{}1,5\] hydrogen shift reaction of all examined compounds can be described as an outcome of two reactions. The second reaction (leading from deuterio-methyl-dideuterio-methylene to dideuterio-methyl-deuterio-methylene in tri-deuterated compounds and from dideuterio-methyl-methylene to deuterio-methyl-deuterio-methylene in di-deuterated compounds) produces a mixture of species that differ by the orientation of the mono-deuterated methylene group. Due to a high barrier for rotation of this group, the two products cannot be properly sampled in a single PIMD simulation. Therefore an additional PIMD simulation is necessary, as shown on the example of **1-5,5,5-*d*$_{\text{3}}$** in Fig. \[fig:PIMD-3steps-pent\]. The reduced free energy of the second step is then calculated as $$\Delta F^{\mathrm{red}}=-k_{B}T\ln\left[\frac{K_{2}^{\mathrm{PIMD}}}{2}\left(1+K_{3}^{\mathrm{PIMD}}\right)\right],\label{eq:F_red_2step_pimd}$$ where 1/2 is the ratio of symmetry factors and $K_{2}^{\mathrm{PIMD}}$ and $K_{3}^{\mathrm{PIMD}}$ stand for equilibrium constants obtained by the second and third PIMD simulations, respectively. Together they represent the second reaction step.

![Three PIMD simulations used to compute equilibrium ratios of the [\[]{}1,5\] hydrogen shift reaction in (3*Z*)-(5,5,5-$^{\text{2}}$H$_{\text{3}}$)penta-1,3-diene (**1-5,5,5-*d*$_{\text{3}}$**). Half white, half black spheres represent deuterium atoms. The methyl group, in contrast to methylene, rotates during simulations. \[fig:PIMD-3steps-pent\]](Fig5_jcp_300dpi)

Combination of *ab initio* methods with the PIMD {#combination-of-ab-initio-methods-with-the-pimd .unnumbered}
------------------------------------------------

Unfortunately, at present the PIMD method cannot be used in conjunction with higher level *ab initio* methods, due to a high number of potential energy evaluations needed. Semiempirical methods, which can be used instead, do not achieve comparable accuracy. We therefore make the following two assumptions: First, we assume that the main contribution to the EIE can be calculated in the framework of HA. Second, we assume that selected semiempirical methods are accurate enough to reliably estimate the anharmonicity correction. The anharmonicity correction is calculated as $$\Delta\Delta F^{\mathrm{anharm}}=\Delta F_{\mathrm{PIMD}}^{\mathrm{red}}-\Delta F_{\mathrm{HA}}^{\mathrm{red}}.\label{eq:F_anharm}$$ With these two assumptions, we can take advantage of both PIMD and higher level methods by adding the semiempirical anharmonicity correction to the HA result calculated by a more accurate method. The HA value of $\Delta F^{\mathrm{red}}$ is obtained by Boltzmann averaging over all possible distinguishable conformations, $${\displaystyle \Delta F_{\mathrm{HA}}^{\mathrm{red}}=-k_{B}T\ln\left[\frac{s^{(r)}\sum_{i=1}^{N_{p}}\left\{ \exp\left(\frac{-E_{i}^{\mathrm{el}}}{k_{B}T}\right)\sum_{j=1}^{s^{(p)}}Q_{ij}^{\mathrm{\mathit{p,\mathrm{nuc}}}}\right\} }{s^{(p)}\sum_{i=1}^{N_{r}}\left\{ \exp\left(\frac{-E_{i}^{\mathrm{el}}}{k_{B}T}\right)\sum_{j=1}^{s^{(r)}}Q_{ij}^{\mathrm{\mathit{r,\mathrm{nuc}}}}\right\} }\right]},\label{eq:F_red_HA}$$ where $N_{r}$ is the number of “geometrically different isomers” of a reactant. By geometrically different isomers we mean species differing in their geometry, not species differing only in positions of isotopically substituted atoms. $E_{i}^{\mathrm{el}}$ is the electronic energy (including nuclear repulsion) of $i^{\mathrm{th}}$ isomer, $s^{(r)}$ is the symmetry factor and $Q_{ij}^{r,\mathrm{nuc}}$ are partition functions of the nuclear motion of $s^{(r)}$ isotopomers. $N_{p}$, $s^{(p)}$, $Q_{ij}^{p,\mathrm{nuc}}$ denote analogous quantities for the product.

Symmetry factors {#symmetry-factors .unnumbered}
----------------

Although symmetry effects can be computed analytically, for the reactions studied in this paper they are nontrivial and so we discuss them here in more detail. As mentioned above, we are interested in moderate temperatures (above $\sim100$ $\mathrm{K}$) where quantum effects might be very important but the distinguishable particle approximation remains valid. In this case, effects of particle indistinguishability and of non-distinguishing several, in principle distinguishable minima, by an experiment can be conveniently unified by the concept of symmetry factor. In our setting, we will call “symmetry factor” the product$$s=s^{\mathrm{exp}}\cdot\prod_{i=1}^{N}\frac{1}{\sigma_{i}}.\label{eq:symmetry_factors}$$ Here, $s^{\mathrm{exp}}$ refers to the number of distinguishable minima not distinguished by the experiment and $\sigma_{i}$ are the usual rotational symmetry numbers ****of symmetric rotors. The symmetry numbers are present only if either the whole molecule or some of its parts are treated as free or hindered classical symmetrical rotors. (In this case, the number of minima of the hindered rotor potential is not included in $s^{\mathrm{exp}}$.) The symmetry numbers are not present if rotational barriers are so high that the corresponding degrees of freedom should be considered as vibrations.

The concept can be illustrated on an example of the mono- and non-deuterated methyl group in a rotational potential with three equivalent minima 120$\,^{\circ}$ apart. At low temperatures, when hindered rotation of the methyl group reduces to a vibration, the symmetry factor is determined only by $s^{\mathrm{exp}}$. In the case of mono-deuterated methyl group there are three in principle distinguishable minima corresponding to three rotamers, which are, as we suppose, considered to be one species by the observer. Therefore $s^{\mathrm{exp}}=3$. In the case of non-deuterated methyl there is only one distinguishable minimum and $s^{\mathrm{exp}}=1$. At higher temperatures, when the methyl group can be treated as a hindered rotor, its contribution to the symmetry factor is determined only by rotational symmetry number $\sigma$, where $\sigma=1$ for a mono-deuterated methyl group and $\sigma=3$ for a non-deuterated methyl group. From definition , it is clear that the high and low temperature pictures are consistent and give the same ratios of symmetry factors.

Partition functions $Q^{(r)}$ and $Q^{(p)}$ needed in calculation of the reduced free energy in Eq. are computed as sums of partition functions of $s^{\mathrm{(r),exp}}$ and $s^{\mathrm{(p),exp}}$ isotopomers.

Computational details\[sec:Computational-details\]
==================================================

Ab initio, density functional, and semiempirical methods {#ab-initio-density-functional-and-semiempirical-methods .unnumbered}
--------------------------------------------------------

In this subsection, we will discuss the accuracy of four electronic structure methods used in our calculations. Ab initio MP2 and the B98 density functional method,[@becke:8554; @schmider:9624] both in combination with the 6-311+(2df,p) basis set, were used for calculations within the HA. Semiempirical AM1[@Dewar1985] and SCC-DFTB[@Elstner1998] methods were used in both HA and PIMD calculations. This allowed us to compute the error introduced by the HA.

Aside from symmetry factors, the EIE is dominated by vibrational contributions. Therefore we concentrate mainly on the accuracy of HA vibrational frequencies. According to Merrick *et al.,*[@Merrick2007] who tested the performance of MP2 and B98 methods by means of comparison with experimental data for a set of 39 molecules, RMSE of ZPEs is 0.46 $\mathrm{kJ\cdot mol^{-1}}$ at MP2/6-311+(2df,p) and 0.31 $\mathrm{kJ\cdot mol^{-1}}$ at B98/6-311+(2df,p) level of theory. Appropriate ZPE scaling factors are equal to 0.9777 and 0.9886 respectively. Corresponding RMSEs of frequencies are 40 $\mathrm{cm^{-1}}$ and 31 $\mathrm{cm^{-1}}.$ Therefore, a slightly higher accuracy can be expected from the B98 functional. The accuracy of both semiempirical methods is significantly worse than the accuracy of higher level *ab initio* methods. According to Witek and Morokuma, the RMSE of AM1 frequencies in comparison to experimental values for 66 molecules is 95 $\mathrm{cm^{-1}}$ with frequency scaling factor equal to 0.9566.[@Witek2004] The error of vibrational frequencies obtained using SCC-DFTB depends on parametrization. We tested two parameter sets: the original SCC-DFTB parametrization[@Elstner1998] and the parameter set optimized with respect to frequencies by Malolepsza *et al.*[@Malolepsza2005] (further designated SCC-DFTB-MWM). The error of vibrational frequencies calculated with the original SCC-DFTB parameters was studied by Krüger *et al.*[@kruger:114110] The mean absolute deviation from the reference values calculated at BLYP/cc-pVTZ for a set of 22 molecules was 75 $\mathrm{cm^{-1}}$. The error of the reference method itself, as compared to an experiment with a slightly smaller set of molecules was 31 $\mathrm{cm^{-1}}$. In the above mentioned study performed by Witek and Morokuma, RMSE of 82 $\mathrm{cm^{-1}}$ with scaling factor 0.9933 was obtained.[@Witek2004] The SCC-DFTB-MWM mean absolute deviation of experimental and calculated frequencies for a set of 14 hydrocarbons is indeed better and is equal to 33 $\mathrm{cm^{-1}}$ instead of 59 $\mathrm{cm^{-1}}$ for the original parametrization. [@Malolepsza2005]

The suitability of AM1, SCC-DFTB, SCC-DFTB-MWM, and several other semiempirical methods for our systems was tested by comparison of EIE values in the HA and of potential energy scans with the corresponding quantities computed with MP2 and B98. Results of this comparison are presented in the appendix. The AM1 method is shown to reproduce *ab initio* EIEs in the HA very well, but fails to reproduce scans of potential energies of methyl and vinyl group rotations. Therefore, in addition to AM1 method we used the SCC-DFTB method, which is, among the semiempirical methods tested by us, the best in reproducing *ab initio* potential energy scans. On the other hand, compared to AM1, SCC-DFTB gives a worse EIE in the HA.

Statistical errors, convergence, and parameters of PIMD simulations {#statistical-errors-convergence-and-parameters-of-pimd-simulations .unnumbered}
-------------------------------------------------------------------

Statistical RMSEs of averages of PIMD simulations were calculated from the equation $$\sigma_{\mathrm{av}}=\tau_{\mathrm{c}}^{1/2}\tau_{\mathrm{sim}}^{-1/2}\sigma,\label{eq:correlated_sigma}$$ where $\sigma$ is the RMSE of one sample. Correlation lengths $\tau_{\mathrm{c}}$ were estimated by the method of block averages.[@flyvbjerg:461] At constant temperature, the correlation length decreases with increasing number of imaginary time slices. Also, for constant number of imaginary time slices, the correlation length decreases with increasing temperature. Finally, the correlation length stays approximately constant at different temperatures, if the number of imaginary time slices is chosen so that $\Delta F$ is converged to approximately the same precision. In our systems, correlation length is close to 3.5 $\mathrm{ps}$.

The EIE was studied at four different temperatures: 200, 441.05, 478.45, and 1000 $\mathrm{K}$ using the normal mode version of PIMD.[@Berne1986] To control temperature, the Nosé-Hoover chains with four thermostats coupled to each PI degree of freedom were used.[@martyna:2635].

Different numbers $P$ of imaginary time slices had to be used at different temperatures, since at lower temperatures quantum effects become more important, and the number of imaginary time slices necessary to maintain the desired accuracy increases. To examine the required value of $P$ as a function of temperature we used the GAFF force field.[@Wang2004] Whereas the accuracy of vibrational frequencies calculated using the GAFF force field is relatively low [\[]{}RMS difference of GAFF and B98/6-311+(2df,p) frequencies of compound **1** is equal to 125 $\mathrm{cm^{-1}}$\], the potential should be realistic enough for the assessment of the convergence with respect to $P$. For example, the difference of potential energies $\Delta E$ between s-*cis* and s-*trans* conformations of compound **1** is 4.3 $\mathrm{kcal\cdot mol^{-1}}$ as compared to 2.7 $\mathrm{kcal\cdot mol^{-1}}$ obtained by MP2/6-311+(2df,p) or 3.5 $\mathrm{kcal\cdot mol^{-1}}$ obtained by B98/6-311+(2df,p).

To check the convergence at 478.45 $\mathrm{K}$, we calculated values of the integral in Eq. for the deuterium transfer reaction in **1-5,5,5-*d*$_{\text{3}}$** with 40 and 48 imaginary time slices (using Simpson’s rule with 5 points). Their difference is equal to 0.00005 $\pm$ 0.00040 $\mathrm{kcal\cdot mol^{-1}}$. This is less than the statistical error of the calculation on the model system, which itself is smaller than the error of production calculations, since the model calculation was 10 times longer than the longest production calculation. Therefore, taking into account the accuracy of production calculations, $P=40$ can be considered the converged number of imaginary time slices. This is further supported by the observation of the convergence of the single value of the GVE at $\lambda=0$. The relative difference of GVE values obtained with 40 and 48 imaginary time slices is equal to 0.22 $\pm$ 0.02 %, whereas the difference between 40 and 72 imaginary time slices equals to 0.39 $\pm$ 0.04 %. The discretization error is asymptotically proportional to $P^{-2}$. By fitting this dependence to the calculated values, we estimated the difference between the values for $P=40$ and for the limit $P\rightarrow\infty$ to be less than 0.6 %. This is only three times more than the difference between $P=40$ and $P=48$. Therefore, if the integral converges similarly as the GVE, we can use the aforementioned difference between 40 and 48 imaginary time slices as the criterion of convergence. The convergence of the GVE with the number of imaginary time slices is displayed in Fig. \[fig:GVE\_TE\_values\]. Since 441.05 $\mathrm{K}$ is close enough to 478.45 $\mathrm{K}$ we used the same value of $P$ at this temperature. To check the convergence at 200 and 1000 $\mathrm{K}$, we observed only the convergence of the single value of the GVE. At 200 $\mathrm{K}$ the relative difference of GVE values at $\lambda=0$ obtained with $P=72$ and $P=80$ is 0.03 $\pm$ 0.07 %. Based on the comparison with the previous result, $P=72$ is considered sufficient. At 1000 $\mathrm{K}$, the relative difference of GVE values at $\lambda=0$ for $P=24$ and $P=32$ equals to 0.1 $\pm$ 0.02 %, so that 24 imaginary time slices are used further.

The time step at 441.05 and 478.45 $\mathrm{K}$ was 0.05 $\mathrm{fs}$ to satisfy the requirement of energy conservation. At 200 and 1000 $\mathrm{K}$, a shorter step of 0.025 $\mathrm{fs}$ was used due to the increased stiffness of the harmonic bonds between beads at 200 $\mathrm{K}$ and due to the increased average kinetic energy at 1000 $\mathrm{K}$. The simulation lengths differed for different molecules. For both isotopologs of compound **1** simulation length of 1 $\mathrm{ns}$ ensured that the system properly explored both the s-*trans* and the s-*cis* conformations. Convergence was checked by monitoring running averages and by comparing the ratio of the s-*trans* and s-*cis* conformers with the ratio calculated in the HA. The length of converged PIMD simulations of compound **2** was 500 $\mathrm{ps}$. Convergence was checked again using running averages and by visual analysis of trajectories to ensure that the system properly explored all local minima. For compound **3,** the simulation length was 400 $\mathrm{ps}$.

The integral in Eq. was calculated using the Simpson’s rule. Using the AM1 potential, the GVE was evaluated for five values of $\lambda$, namely for $\lambda=0.0,0.25,0.5,0.75,1.0$. Convergence was checked by comparison with values obtained using the trapezoidal rule. Since the dependence of the estimator on the parameter $\lambda$ is almost linear, the difference between the two results remained well under the statistical error. Using the SCC-DFTB potential, the dependence on $\lambda$ was less smooth. As a result, nine equidistant values of $\lambda$ were needed to achieve similar convergence. In one case, as many as 17 values of $\lambda$ had to be used.

[lcccc]{} & 1$^{\text{st}}$ step $\Delta F^{\mathrm{red}}$ & $\Delta\Delta F^{\mathrm{anharm}}$ & 2$^{\text{nd}}$ step $\Delta F^{\mathrm{red}}$ & $\Delta\Delta F^{\mathrm{anharm}}$ [\
]{} [\
]{} AM1 (PIMD)& 0.0395& -0.0041$\pm$0.0009& -0.0154& 0.0022$\pm$0.0007[\
]{}SCC-DFTB (PIMD)& 0.1245& -0.0039$\pm$0.0007& -0.0616& 0.0026$\pm$0.0005[\
]{}B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$& 0.0587& & -0.0248& [\
]{}MP2 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$& 0.0770& & -0.0338& [\
]{} [\
]{}AM1 (PIMD)& -0.0283& 0.0063$\pm$0.0009& 0.0191& -0.0023$\pm$0.0007[\
]{}SCC-DFTB (PIMD)& -0.1142& 0.0049$\pm$0.0006& 0.0610& -0.0017$\pm$0.0005[\
]{}B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$& -0.0466& & 0.0282& [\
]{}MP2 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$& -0.0645& & 0.0372& [\
]{}

Amber 10 implementation {#amber-10-implementation .unnumbered}
-----------------------

The PIMD calculations were performed using Amber 10.[@amber10] The part of Amber 10 code, which computes the derivative $dF\left(\lambda\right)/d\lambda$ with respect to the mass was implemented by one of us and can be invoked by setting ITIMASS variable in the input file. Several posible ways to compute the derivative are obtained by combining one of two implementations of PIMD in `sander (`either the multisander implementation or the LES implementation) with either the Nosé-Hoover chains of thermostats or the Langevin thermostat, with the normal mode or “primitive” PIMD, and with the TE or GVE. Calculating the value of $dF\left(\lambda\right)/d\lambda$ for compound **1-5,5,5-*d*$_{\text{3}}$** using GAFF force field, at $\lambda=0$, $T=478.45$ $\mathrm{K}$, and for several values of $P$, we confirmed that all twelve possible combinations give the same result. For example, Fig. \[fig:GVE\_TE\_values\] shows the agreement of the GVE and TE. Nevertheless, the twelve methods differ by RMSEs of $dF\left(\lambda\right)/d\lambda$ (due to different statistical errors of estimators and correlation lengths) and by computational costs. As expected, the most important at higher values of $P$ is the difference of statistical errors of GVE and TE. Figure \[fig:GVE\_TE\_errors\] compares the dependence of the RMSE of the GVE and TE on the number $P$ of imaginary time slices. For $P=64$, the converged simulation with GVE is approximately 100 times faster than with TE. Less significant differences of RMSEs are due to differences in correlation lengths. As expected, for primitive PIMD with Langevin thermostat the correlation length depends strongly on collision frequency $\gamma$ of the thermostat. The correlation length is approximately 450-500 $\mathrm{fs}$ for $\gamma=0.3$ $\mathrm{ps^{-1}}$, falling down quickly to 120-150 $\mathrm{fs}$ for $\gamma=3$ $\mathrm{ps^{-1}}$ and to 5 $\mathrm{fs}$ for $\gamma=300$ $\mathrm{ps^{-1}}$, then rising slowly again. This can be compared to correlation length of 10-20 $\mathrm{fs}$ of the primitive PIMD thermostated by Nosé-Hoover chains of 4 thermostats per degree of freedom. A smaller difference can be found between correlation lengths of the normal mode (3-8 $\mathrm{fs}$) and primitive PIMD (10-20 $\mathrm{fs}$).

Used software {#used-software .unnumbered}
-------------

All PIMD calculations were performed in Amber 10.[@amber10] All MP2 and B98 calculations as well as AM1 and PM3 semiempirical calculations in the HA were done in Gaussian 03 revision E01.[@g03] DFTB and SCC-DFTB calculations in the HA used the DFTB+ code, version 1.0.1.[@AradiB._jp070186p] SCC-DFTB harmonic frequencies were computed numerically using analytical gradients provided by the DFTB+ code. The step size for numerical differentiation was set equal to 0.01 $\textrm{\AA}$. This value was also used by Krüger *et al.*[@kruger:114110] in their study validating the SCC-DFTB method and their frequencies differed from purely analytical frequencies of Witek *et al.*[@Witek2004] by at most 10 $\mathrm{cm^{-1}}$. To diagonalize the resulting numerical Hessian, we used the `formchk` utility included in the Gaussian program package.

Results\[sec:Results\]
======================

(3*Z*)-penta-1,3-diene
----------------------

(3*Z*)-penta-1,3-diene (compound **1**) is the simplest of examined molecules. Its non-deuterated isotopolog has three distinguishable minima: the s-*trans* conformer, which is the global minimum and two s-*cis* conformers related by mirror symmetry. Strictly speaking, in pentadiene the s-*cis* species have actually *gauche* conformations due to sterical constraints. In their original experiments, Roth and König studied two isotopologs, **1-5,5,5-*d*$_{\text{3}}$** and **1-1,1-*d*$_{\text{2}}$**. The EIEs of both isotopologs were computed using the PIMD methodology of Sec. \[sec:Method\] at 478.45 $\mathrm{K}$. Resulting reduced free energies are listed in Table \[tab:F\_pentadiene\]. Note that the anharmonicity correction is very similar for AM1 and SCC-DFTB methods, even though the main part of obtained in the HA ( substantially differs for the two methods. This indicates that the corrections are fairly reliable in this case and can be used to correct results of higher level methods obtained in the HA. The anharmonicity correction is as large as 20 % of the final value of the reduced free energy. Unfortunately, the difference between MP2 and B98 in the HA is still approximately four times larger than the anharmonicity correction.

                                                                  **1-5,5,5-*d*$_{\text{3}}$**   **1-1,1,5-*d*$_{\text{3}}$**   **1-1,5,5-*d*$_{\text{3}}$**
  -------------------------------------------------------------- ------------------------------ ------------------------------ ------------------------------
  AM1 (PIMD)                                                                 0.103                          0.296                          0.601
  SCC-DFTB (PIMD)                                                            0.108                          0.285                           0606
  B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$               0.104                          0.293                          0.602
  MP2 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$               0.105                          0.291                          0.604
                                                                   **1-1,1-*d*$_{\text{2}}$**     **1-5,5-*d*$_{\text{2}}$**     **1-1,5-*d*$_{\text{2}}$**
  AM1 (PIMD)                                                                 0.099                          0.305                          0.597
  SCC-DFTB (PIMD)                                                            0.093                          0.315                          0.591
  B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$               0.097                          0.307                          0.596
  MP2 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$               0.096                          0.309                          0.595

  : Equilibrium ratios of [\[]{}1,5\] hydrogen shift reactions of **1-5,5,5-*d*$_{\text{3}}$** and **1-1,1-*d*$_{\text{2}}$** at 478.45 $\mathrm{K}$. \[tab:eq\_ratios\_pentadiene\]

                                                                  1$^{\text{st}}$ step $\Delta F^{\mathrm{red}}$   $\Delta\Delta F^{\mathrm{anharm}}$   2$^{\text{nd}}$ step $\Delta F^{\mathrm{red}}$   $\Delta\Delta F^{\mathrm{anharm}}$
  -------------------------------------------------------------- ------------------------------------------------ ------------------------------------ ------------------------------------------------ ------------------------------------
  AM1 (PIMD)                                                                         -0.0337                               0.0102$\pm$0.0013                                0.0230                               -0.0036$\pm$0.0010
  B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$                       -0.0496                                                                                0.0302                      
  MP2 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$                       -0.0674                                                                                0.0392                      

Reduced free energies suggest a general preference valid for both tri- and di-deuterio species: Namely, deuterium, compared to hydrogen, prefers sp$^{\text{3}}$ carbon of the methyl group to sp$^{\text{2}}$ carbon of the vinyl group. This was also observed experimentally.[@SunkoDionisE._ja00725a026; @Barborak1971; @Gajewski1979] As can be seen in the table, the preference is present already in the HA. At first sight this preference can be counter-intuitive: If we confine ourselves only to the most energetical C-H (or C-D) bond stretching modes and suppose that force constants do not change significantly upon substitution with deuterium, the deuterium should prefer the stiffest bonds. This is because more energy can be gained by substituting the stiffer sp$^{2}$ C-H bonds, assuming approximately the same change in reduced mass after substitution. (Recall that the energy of a vibrational mode is proportional to $\sqrt{k/\mu}$ where $k$ is the force constant and $\mu$ the reduced mass.) Considering the stretching modes only, this would be the case for isotopologs of compound **1**. Nevertheless, taking into account also the bending and torsional vibrational modes, gains on the sp$^{\text{3}}$ C-H bond side will dominate. (This “counter-intuitive” preference of heavier isotope in “softer” bonds is quite common. For examples see the inverse H/D equilibrium isotope effect in oxidative addition reactions of H$_{\text{2}}$ to transition metal complexes,[@Hascall1999; @Bender1995; @Abu-Hasanayn1993; @Rabinovich1993] or the inverse $^{\text{16}}$O/$^{\text{18}}$O isotope effect in metal mediated oxygen activation reaction.[@Smirnov2009]) As already stressed, final equilibrium concentrations are determined mainly by the symmetry factors and the aforementioned deuterium sp$^{\text{3}}$ to sp$^{\text{2}}$ preference manifests itself only in a small modification of the symmetry determined rational ratios, as seen in Table \[tab:eq\_ratios\_pentadiene\].

### Temperature dependence of the reduced free energy {#temperature-dependence-of-the-reduced-free-energy .unnumbered}

Temperature dependence of the reduced free energy for the first reaction step of hydrogen shift in **1-5,5,5-*d*$_{\text{3}}$** is depicted in Fig. \[fig:harmonic\_freq\_semiempirical\_test\]. Analogous temperature dependence for all other studied reactions is very similar to the dependence in this figure. At very low temperature, the reduced reaction free energy approaches the difference of ZPEs in accordance with Eq. . Absolute value of the reduced free energy is maximal at temperatures around 200 $\mathrm{K}$. At temperatures around 400-500 $\mathrm{K}$, where most measurements took place, the value of is (by chance) close to the difference of ZPEs. At high temperatures, goes to zero in accordance with the high temperature limit discussed above, which is valid also in the HA as can be shown using the Teller-Redlich theorem.[@Angus1935; @Redlich1935] The temperature dependence of the anharmonicity correction was examined using the AM1 semiempirical method at temperatures 200 $\mathrm{K}$, 478.45 $\mathrm{K}$ and 1000 $\mathrm{K}$. The correction for the first reaction step of **1-5,5,5-*d*$_{\text{3}}$** is equal to $-0.0042$ $\pm$ $0.0009$, $-0.0041$ $\pm$ $0.0009$ and $-0.0035$ $\pm$ $0.0008$ $\mathrm{kcal\cdot mol^{-1}}$, respectively. Therefore, taking into account statistical errors, the correction stays approximately constant over the wide temperature range. Since the value of is decreasing in the region from 200 to 1000 $\mathrm{K}$, the relative importance of the correction is increasing.

![Reduced reaction free energies of the first step of [\[]{}1,5\] hydrogen shift reaction in (3*Z*)-(5,5,5-$^{\text{2}}$H$_{\text{3}}$)penta-1,3-diene (**1-5,5,5-*d*$_{\text{3}}$**) calculated in the HA as the Boltzmann average of all s-*trans* and s-*cis* isomers.\[fig:harmonic\_freq\_semiempirical\_test\]](Fig6_jcp_300dpi)

2-methyl-10-methylenebicyclo[\[]{}4.4.0\]dec-1-ene 
---------------------------------------------------

                                                                  **2-1,1-*d*$_{\text{2}}$**   **2-5,5-*d*$_{\text{2}}$**   **2-1,5-*d*$_{\text{2}}$**
  -------------------------------------------------------------- ---------------------------- ---------------------------- ----------------------------
  AM1 (PIMD)                                                                0.098                        0.306                        0.595
  B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$              0.097                        0.308                        0.595
  MP2 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$              0.096                        0.310                        0.594
  experimental (series 1)                                                   0.108                        0.328                        0.564
  experimental (series 2)                                                   0.114                        0.314                        0.572

  : Equilibrium ratio of the [\[]{}1,5\] hydrogen shift reaction in dideuterated compound **2** at 441.05 $\mathrm{K}$. Experimental series 1 and 2 were obtained by two different methods of analysis of the NMR spectrum.[@Doering2006][]{data-label="tab:eq_ratios_bicyclo"}

                1$^{\text{st}}$ step $\Delta F^{\mathrm{red}}$   $\Delta\Delta F^{\mathrm{anharm}}$   2$^{\text{nd}}$ step $\Delta F^{\mathrm{red}}$   $\Delta\Delta F^{\mathrm{anharm}}$
  ------------ ------------------------------------------------ ------------------------------------ ------------------------------------------------ ------------------------------------
  AM1 (PIMD)                       -0.0439                               0.0160$\pm$0.0015                                0.0265                               -0.0080$\pm$0.0011
                                   -0.0673                                                                                0.0376                      

Compound **2-1,1-*d*$_{\text{2}}$** was used relatively recently by Doering and Zhao[@Doering2006] to refine the original results of Roth and König. Doering and Zhao reported the equilibrium concentrations at three temperatures, from which we have chosen the lowest, $T=441.05$ $\mathrm{K}$. Since AM1 and SCC-DFTB methods gave similar anharmonicity corrections for compound **1**, we used only the AM1 method in this case. Four minima were found. Energy differences between the global minimum and local minima at the B98/6-311+(2df,p) level of theory are 4.8, 7.5 and 9.5 $\mathrm{kcal\cdot mol^{-1}}$. Resulting calculated according to Eq. are listed in Table \[tab:F\_bicyclo\]. As expected, they are similar to of **1-1,1-*d*$_{\text{2}}$**. Anharmonicity corrections changed more and they are about 60 % higher in the absolute value. Calculated equilibrium ratios are listed in Table \[tab:eq\_ratios\_bicyclo\] together with experimental ratios reported by Doering and Zhao. Theoretical and experimental ratios differ substantially, which suggests that side reactions suspected by Doering and Zhao had indeed occurred and influenced the accuracy of results of their study.

2,4,6,7,9-pentamethyl-5-methylene-11,11a-dihydro-12*H*-naphthacene
------------------------------------------------------------------

In order to suppress unwanted side reactions suspected for compound **2**,[@Doering2006] Doering and Keliher further developed the model compound into **3-1,1-*d*$_{\text{2}}$** by adding methyl substituted aromatic rings on both sides of the cyclic part.[@Doering2007] With this compound they obtained the same equilibrium ratios for all temperatures they had examined.[@Doering2007] Because of this and because the temperature 441.05 $\mathrm{K}$ we have chosen for our analysis of compound **2** differs from one of temperatures used in Ref. by less than 2 $\mathrm{K}$, we have decided to use this temperature also for compound **3**. In the HA, only the B98 density functional method was used, due to a considerable size of the molecule. Because of the increased rigidity imposed by aromatic rings on the sides of the original bi-cyclic compound, only three distinct minima were found. (We neglected several possible orientations of two methyl groups distant from the reaction site, which hardly affect the final result.) At the B98/6-311+G(2df,p) level, the local minima have energies 1.3 and 5.6 $\mathrm{kcal\cdot mol^{-1}}$ above the global minimum. The AM1 method gives the opposite order of the first and second lowest minima. Reduced free energies and anharmonicity corrections obtained using the AM1 method are listed in Table \[tab:F\_naphthacene\].

Values of both and the anharmonicity corrections are again qualitatively similar (but higher in absolute values) to those of the smaller and less strained compound **2**. Here, values of anharmonicity corrections reached approximately 30 % of the values of reduced free energies of both reaction steps. Resulting equilibrium ratios together with their experimental values can be seen in Table \[tab:eq\_ratios\_naphthacene\]. Agreement of the theoretical prediction with the experimental result is very good. An uncorrected B98 HA value is also included in Table \[tab:eq\_ratios\_naphthacene\], to demonstrate how the anharmonicity correction modifies the HA equilibrium ratio. Surprisingly, the direct AM1 PIMD calculation is closest to the experimental value, but this ought to be ascribed to a fortunate coincidence, considering the aforementioned accuracy of electronic structure methods used in our study.

                                                                  **3-1,1-*d*$_{\text{2}}$**   **3-5,5-*d*$_{\text{2}}$**   **3-1,5-*d*$_{\text{2}}$**
  -------------------------------------------------------------- ---------------------------- ---------------------------- ----------------------------
  AM1 (PIMD)                                                                0.098                        0.308                        0.595
  B98 (HA)                                                                  0.095                        0.312                        0.593
  B98 (HA) + $\Delta\Delta F_{\mathrm{AM1}}^{\mathrm{anharm}}$              0.096                        0.310                        0.594
  experimental                                                              0.098                        0.308                        0.594

  : Equilibrium ratios of the [\[]{}1,5\] hydrogen shift reaction in compound **3** at 441.05 $\mathrm{K}$. []{data-label="tab:eq_ratios_naphthacene"}

Discussion and conclusions\[sec:Discussion-and-conclusions\]
============================================================

To conclude, the combination of higher level methods in the HA with PIMD using semiempirical methods for the rigorous treatment of effects beyond the HA proved to be a viable method for accurate calculations of EIEs. Using the generalized virial estimator for the derivative of the free energy with respect to the mass we were able to obtain accurate results at lower temperatures in reasonable time ($\sim60$ times faster than with thermodynamic estimator), since the statistical error is independent on the number of imaginary time slices. Two semiempirical methods, AM1 and SCC-DFTB, were used for calculation of the anharmonicity correction, both giving very similar results. Calculations showed, that the anharmonicity effects account up to 30 % of the final value of the reduced free energy of considered reactions. The anharmonicity correction always decreases the absolute value of the reduced reaction free energy. This is consistent with the qualitative picture in which higher vibrational frequencies of hydrogens are lowered by the anharmonicity of a potential more than lower frequencies of deuteriums. This in turn is due to a higher amplitude of vibrations of lighter hydrogen atoms. The lower difference between frequencies of unsubstituted and deuterated species results in the lower absolute value of the reduced reaction free energy. Unfortunately, the inaccuracy of the *ab initio* electronic structure methods used in our study is still of at least the same order as the anharmonicity corrections.

For isotopologs of compound **1**, we predicted equilibrium ratios and free energies of the [\[]{}1,5\] sigmatropic hydrogen shift reaction. A comparison with experimental results was not possible due to a low precision of the original measurement. For compound **2**, the disagreement between theoretical and experimental data supports the suspicion by authors of the measurement that the accuracy of their results was compromised by dimerization side reactions. On the other hand, the agreement of theoretically calculated ratios with experimental observations in the case of compound **3** suggests that the isolation of the [\[]{}1,5\] hydrogen shift reaction from disturbing influences was successfully achieved and the observed EIE and KIE can be considered reliable.

We acknowledge the start-up funding provided by EPFL. We express our gratitude to Daniel Jana for his assistance with SCC-DFTB calculations, to H. Witek for providing us the set of frequency optimized SCC-DFTB parameters, and to J. J. P. Stewart for the MOPAC2007 code used for testing the PM6 method.

Examination of semiempirical methods used in PIMD simulations
=============================================================

To determine the suitability of semiempirical methods used in our PIMD calculations, we first compared values of the EIE in the HA. The B98 and MP2 methods served as a reference. As can be seen from Fig. 6 in the paper, which shows the temperature dependence of $\Delta F^{\mathrm{red}}$ of the first step of deuterium transfer reaction in **1-5,5,5-*d*$_{\text{3}}$**, the difference between B98 and MP2 at temperatures below 500 $\mathrm{K}$ is close to 0.02 kcal/mol. So is the difference between AM1 and B98 methods. On the other hand, the SCC-DFTB method clearly overestimates the extent of EIE compared to both higher level methods. A very similar trend was observed in all examined reactions. Other semiempirical methods tested were PM3[@Stewart1989] and SCC-DFTB-MWM, which are not included in Fig. 6 in the paper for clarity. The PM3 method overestimates the EIE similarly to SCC-DFTB, whereas the SCC-DFTB-MWM curve is somewhat closer to *ab initio* curves than the SCC-DFTB one. To conclude, from this point of view AM1 is the preferred semiempirical method.

During simulations, the pentadiene molecule often passes two potential energy barriers. These are the barrier for the hindered rotation of the methyl group and the barrier for the rotation of the vinyl group, which connects s-*trans* and s-*cis* conformations. Relaxed potential energy scans of these two motions were employed as the second criterion to assess the relevancy of semiempirical methods. Methods tested were MP2, B98, AM1, SCC-DFTB, SCC-DFTB-MWM, PM3, RM1,[@Rocha2006] PM3CARB-1,[@McNamara2004] PDDG/PM3,[@Repasky2002] and PM6.[@Stewart2007] Potential surface scans with PM3CARB1, RM1, PM3/PDDG methods were calculated using the public domain code MOPAC 6. PM6 potential surface scans were performed in MOPAC 2007.[@mopac2007] Results for the methyl group rotation are shown in Fig. \[fig:methyl\_scan\]. The height of the AM1 barrier is only 0.005 $\mathrm{kcal\cdot mol^{-1}}$. Moreover, positions of minima do not agree with B98 and MP2. On the other hand, the SCC-DFTB method matches higher level methods closely. From other semiempirical methods PM3 performs best in this aspect. The height of the barrier is relatively well reproduced also by PDDG/PM3, RM1 and SCC-DFTB-MWM methods, but positions of extrema of the potential energy surface are incorrect.

![Relaxed potential energy scan of the methyl rotation in s-*trans* (3*Z*)-penta-1,3-diene (compound **1**). For MP2, B98, DFTB, SCC-DFTB and SCC-DFTB-MWM only positions and potential energies of minima and maxima are indicated. \[fig:methyl\_scan\]](Fig_supp_info_1_jcp_300dpi)

Figure \[fig:vinyl\_scan\] shows potential energy scans of the s-*trans*/s-*cis* rotation of the vinyl group. Again, the SCC-DFTB method matches higher level methods closely. All other methods (with the exception of DFTB) give too low barrier heights as well as too low energy differences between s-*trans* and s-*cis* conformations. Also note that the potential energy surfaces of PM3 and related methods (PDDG/PM3 and PM3CARB-1) are not smooth in the gauche region. This peculiarity of the PM3 potential surface can be seen also in the potential surface scan performed by Liu *et al.*[@Liu1993]

![Relaxed potential energy scan of the vinyl group rotation in s-*trans* (3*Z*)-penta-1,3-diene (compound **1**). For MP2, B98, DFTB, SCC-DFTB and SCC-DFTB-MWM only positions and potential energies of minima and transition states are indicated.\[fig:vinyl\_scan\]](Fig_supp_info_2_jcp_300dpi)

Based on these results we concluded that none of the semiempirical methods except for SCC-DFTB is able to sufficiently improve the AM1 potential energy surface. Whereas the frequency optimized variant of the SCC-DFTB method (SCC-DFTB-MWM) improves the EIE in the HA, it does not retain the SCC-DFTB accuracy in the potential surface scans. Hence we decided to use the AM1 and SCC-DFTB potentials for PIMD calculations. To conclude, AM1 performs very well in HA, but it cannot properly describe potential surfaces of the two rotational motions realized during simulations. On the other hand, SCC-DFTB gives worse results in HA, but it reproduces both barriers very well.

[68]{} natexlab\#1[\#1]{}bibnamefont \#1[\#1]{}bibfnamefont \#1[\#1]{}citenamefont \#1[\#1]{}url \#1[`#1`]{}urlprefix\[2\][\#2]{} \[2\]\[\][[\#2](#2)]{}

, , , , , ****, ().

, , , ****, ().

, , , ****, ().

, ****, ().

, , , , ****, ().

, ****, ().

, , , ****, ().

, , , , , ****, ().

, ****, ().

, , , , ****, ().

, , , , ****, ().

, , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, **, Los Alamos Series in Basic and Applied Sciences (, ).

, ****, ().

, ****, ().

, ** (, , ).

, **, International Series in Pure and Applied Physics (, ).

, , , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , ****, ().

, in **, edited by , , (, ).

, ****, ().

, , , , , , , , , , , ** (), .

, , , , , , , , , , , ** (), .

, ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , , , , , , ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, , , , ****, ().

, ****, ().

, ****, ().

, , , ****, ().

, , , , , ****, ().

, , , , , , , , , , , **, .

, , , ****, ().

, , , , ****, ().

, , , ****, ().

, ****, ().

, , , , , , ****, ().

, ****, ().

, , , ****, ().

, ****, ().

, , , , , , , , ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , , , , , ****, ().

, , , ****, ().

, ****, ().

, ** (), .

, , , , , , ****, ().

[^1]: Namely, the barrier height of internal hindered rotation of the methyl group.
