FLASH INFORMATIQUE FI

High Performance Computing


Path integral evaluation of the equilibrium isotope eff ect in [1,5] sigmatropic hydrogen shift reactions




Tomáš ZIMMERMANN

Jiri VANÍČEK


The path integral methodology for computation of the equilibrium isotope effect is described and applied to three [1,5] sigmatropic hydrogen shift reactions. An efficient estimator for the derivative of the free energy is used, which shortens the computation time by the factor of about 60. Calculations reveal that the correction beyond the usual harmonic approximation for nuclear motion 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 as well as concerns about compromised accuracy, due to side reactions, of another measurement.

[1,5] sigmatropic hydrogen shift reaction

[1,5] sigmatropic hydrogen shift reaction is a molecular rearangment that involves both the creation of a new σ-bond between hydrogen atom and the vinyl carbon and the breaking of an existing π-bond between hydrogen and methyl carbon atom accompanied by a concurrent relocation of π-bonds (see Figure 1).

JPEG - 7.5 kb
fig.1
The [1,5] sigmatropic hydrogen shift reaction in (3Z)-penta-1,3-diene. The reactant and product are identical for the isotopically unsubstituted compound but they can differ when some hydrogens are substituted by deuteriums.

The barrier to this reaction is high so the experimental measurements are usually performed at temperatures around 160-200°C. The reaction attracts attention of both experimental and theoretical chemists for reason that the quantum tunneling supposedly contributes significantly to the reaction rate. Therefore, simple compounds capable of [1,5] hydrogen shift can serve as models to understand phenomena important, for example, in much more complex biological systems.
One approach to answer the question whether the tunneling is significant in the reaction is to study the kinetic isotope effect (KIE). The primary H/D KIE describes how the reaction rate changes when the hydrogen atom being transferred is substituted by deuterium. A high value of the KIE then means that the quantum tunneling can be an important pathway in the reaction. Usually, not only the hydrogen being transferred but all hydrogens bonded to the same carbon atom are replaced by deuteriums. If a molecule studied is symmetric with respect to the [1,5] hydrogen shift (as (3Z)-penta-1,3-diene shown in Figure 1), the reaction products differ only in positions of deuteriums. Their equilibrium ratio is the special case of the equilibrium (thermodynamic) isotope effect (EIE).
Recently, two subsequent measurements of the KIE of [1,5] sigmatropic hydrogen shift were published [1, 2]. For both measurements, the final equilibrium ratios of isotopomers were reported, which were, despite of the similarity of compounds used, qualitatively differed. The authors of the measurements suspected that the hydrogen shift reaction in the first measurement was influenced by the dimerization side reactions.

Equilibrium isotope effect

The EIE is defined as the effect of isotopic substitution on the equilibrium constant. Denoting an isotopolog [1] with a lighter (heavier) isotope by a subscript l (h), the EIE is defined as the ratio of equilibrium constants
EIE = \frac{K^i} {K^h} = \frac{{Q_{l}^{(p)}}/{Q_{l}^{(r)}}}{{Q_{h}^{(p)}}/{Q_{h}^{(r)}}},       (1)
where Q(r) and Q(p) are molecular partition functions of the reactant and product, respectively. We study a specific case of the EIE - the equilibrium ratio of two isotopomers [2]. In this case, the EIE is equal to the equilibrium constant of the isotopomerization reaction,
EIE = K_eq = \frac{Q^{(p)}}{Q^{(r)}},       (2)
where superscripts r and p refer to the reactant and product isotopomers, respectively.

Thermodynamic integration

The direct computation of the ratio of partition functions in Eq. (2) is, in general, difficult, for the value of the partition function is proportional to the huge number of energy levels accessible to the system at a given temperature. To make the computation feasible, a trick of the thermodynamic integration with respect to the mass of the isotopes is used. The method takes advantage of the relationship
EIE = \frac {Q^{(p)}}{Q^{(r)}} = \int_0^{1} d\lambda \frac {dF(\lambda)}{d(\lambda)},       (3)
where F=-logQ/β is the (quantum) free energy and λis a parameter which provides a smooth transition between isotopomers r and p.
In contrast to the partition function itself, the integrand dF(\lambda)}{d{\lambda}} is a thermodynamic average and therefore can be computed by either Monte Carlo or molecular dynamics simulations. We use the PIMD method, which combines the path integral expression for the partition function with molecular dynamics for sampling of the configuration space.

Path integral method

Since the electronic properties of hydrogen and deuterium are the same (to a high level of approximation), the main difference between isotopomers is due to the motion of the nuclei, which move in the same potential but differ in their masses. To go beyond usual approximations and describe the quantum effects on the nuclear motion rigorously, we use the path integral (PI) formulation of quantum mechanics [3]. In the path integral formalism, thermodynamic properties are computed exploiting the correspondence between matrix elements of the Boltzmann operator and the quantum propagator in imaginary time. The PI representation expression for the partition function Q is
Q = Cdr(0)dr(P-1) exp[- β Φ (r(s))],        (4)
where
C = \left( \frac{P}{2\rho h^2 b} \right)^{NPD/2}   \prod\limits_{i=1}^N {m_{i}{pd/2}},        (5)
and  {\bold{r}^(s)} = ({{\bold{r}_1^{(s)}},\bold{r}_2^{(s)}, ... \bold{r}_N^{(s)}}) is the set of Cartesian coordinates associated with the sth imaginary time slice,

\Phi {\left( \left{ \bold{r}^s \right} \right)} = {\frac {P}{2{h^2}{b^2}}} {\sum_{i=1}^{N}} m_i {\sum_{s=0}^{p-1}} \left( \bold{r}_i ^{(s)} - \bold{r}_i ^{(s+1)}\right)^2  + {\frac{1}{P}}{\sum_{s=0}^{p-1}} V \left( \bold{r}^{(s)} \right),    (6)
is the effective potential, N is 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).
The P particles representing each nucleus in P different imaginary time slices are called beads. Each bead interacts with the two beads representing the same nucleus in adjacent time slices via the harmonic potential and with beads representing other nuclei in the same imaginary time slice via the potential of the molecule attenuated by the factor 1/P. See Figure 2 showing the classical system corresponding to a quantum mechanical description of equilibrium of (3Z)-penta-1,3-diene in the case of P = 4.

JPEG - 5.5 kb
fig. 2
The path integral representation of the quantum partition function of (3Z)-penta-1,3-diene for P=4 imaginary time slices represented by different colors. Harmonic interaction between beads in adjacent time slices is represented by the springs. In addition to this harmonic interaction, beads of the same color interact by the potential of the molecule diminished by the factor 1/P. During the simulation, the whole system evolves according to the classical dynamics

The number of imaginary slices, which have to be used to obtain converged quantum result, depends on the system and the temperature. At higher temperature, a system generally behaves more classically, so a lower number P of imaginary time slices is needed. Since the EIE is a small quantity, we mainly used relatively high number of P = 40 imaginary time slices to obtain precise results.

Thermodynamic and generalized virial estimators

To use the PI expression for the partition function together with the thermodynamic integration, we have to calculate the logarithmic derivative [recall Eq. 3] of the partition function Q with respect to the mass. By straightforward differentiation of Eq. (4) , we obtain the so-called thermodynamic estimator (TE),

\frac{dF(\lambda)}{d{\lambda}}
 \simeq 
\frac{-1}{b} \sum_{i=1}^{N} \frac{dm_i}{d\lambda} 
\left< {
{\frac{D,P}{2m_i}} - {\frac{P}{2h^2b}}  \sum_{i=0}^{P-1}
\left(\bold{r}_i ^{(s) - \bold{r}_i ^{(s+1)}\right)^2
 }\right>,       (7)
The problem of the TE is that its statistical error grows with the number P of time slices of the PI, therefore much longer simulation is needed, if one wants to calculate precise result close to the quantum limit.
Fortunately, the growth of the statistical error can be removed if the centroid coordinate of the beads representing the same atom is subtracted, and the coordinates are mass-scaled, prior to performing the derivative [4]. The error of the resulting generalized virial estimator (GVE) stays approximately independent of the number of imaginary time slices
\frac{dF(\lambda)}{d{\lambda}}
 \simeq 
\frac{-1}{\beta} \sum_{i=1}^{N} \frac{dm_i d\lambda}{m_i} 
 x 
\left[ {\frac{D}{2} - \frac{\beta}{2P} \left< \sum_{i=0}^{P-1}\left(\bold{r}_i ^{(s) - \bold{r}_i ^{(C)}\right) \frac{\theta V (\bold {r}^{(s)}}{\theta \bold {r}_i^{(s)} }\right>
}\right] .       (8)
The dependence of the statistical error of both estimators on the number of imaginary time slices in the PI is shown in Figure 3. In the systems studied in this work, the GVE shortened the computations by a factor of about 60 compared to the TE.

JPEG - 6.9 kb
fig. 3
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.

Combination of PIMD with single point ab initio methods

Unfortunately, nowadays the PIMD method cannot be used directly in conjunction with highly precise ab initio methods for calculation of the potential energy surface, due to a high number of energy evaluations needed. Semiempirical methods, which can be used instead, generally do not achieve comparable accuracy. We therefore made the following two assumptions:
First, we assume that the main contribution to the EIE can be calculated by a single point ab initio calculation in the framework of the harmonic approximation (HA). In the HA, the multidimensional potential surface on which nuclei move is approximated by a parabolic potential centered at the local minimum, so the vibrational motion of nuclei is described by the system of coupled harmonic oscillators for which the expression for the EIE [Eq. (2)] can be obtained analytically.
Second, we assume that carefully selected semiempirical methods are accurate enough to reliably estimate the error due to the anharmonicity of the real potential.
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.

Computational details

All computations were done on the cluster of the LCPT group, which consists of 32 nodes each with two Quad-Core AMD Opteron 2352 processors.
The PIMD calculations were performed using Amber 10 [5]. The part of the Amber 10 code, which computes the derivative with respect to the mass was implemented by one of us and can be invoked by setting the ITIMASS variable in the input file. PIMD simulations are well parallelized; the dynamics of each imaginary time slice runs on a single core. Still, the computation of one value of the logarithmic derivative of the partition function [see Eq. (3)] needed to compute the equilibrium ratio of two isotopomers of the largest molecule studied (2,4,6,7,9-pentamethyl-5-methylene-11,11a-dihydro-12H-naphthacene) with the number of imaginary time slices P = 40 took approximately 9 days. Five values of the derivative were usually sufficient to calculate the integral in Eq. (3).
Ab initio single point calculations were done using Gaussian 03 revision E01 [6]. Due to the limitations of the Gaussian 03 parallelization, all ab initio computations were parallelized only on a single node using SMP. Nevertheless, finding local minima and harmonic frequencies for the largest studied molecule lasted only five days instead of one month.

Results

The reaction free energies and equilibrium ratios for two isotopologs of the deuterated (3Z)-penta-1,3-diene and two related compounds were computed. Here, we shortly summarize main results for the latter two compounds. More details and results for (3Z)-penta-1,3-diene can be found in Ref. [6].

JPEG - 6 kb
fig. 4
Global minimum of 2-methyl-10-methylenebicyclo[4.4.0]dec-1-ene

2-methyl-10-methylenebicyclo[4.4.0]dec-1-ene

JPEG - 8.5 kb
fig. 5
Global minimum of 2,4,6,7,9-pentamethyl-5-methylene-11,11a-dihydro-12H-naphthacene

Theoretical and experimental [7] ratios for 2-methyl-10-methylenebicyclo[4.4.0]dec-1-ene (Figure 4) differ substantially, which suggests that side-reactions suspected by Doering and Zhao had indeed occurred and compromised the accuracy of the measurement of the kinetic isotope effect.

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

Theoretical and experimental ratios of products of the [1,5] sigmatropic hydrogen shift in di-deuterated compound (see Figure 5) agree very well as can be seen in Table 1. This agreement supports the accuracy of the latter measurement of the kinetic isotope effect (KIE) made on this compound.

<7h>
5,5-d26,6-d25,6-d2
computed0.0960.3100.594
experimental (Ref. [1])0.0980.3080.594

Table 1 - Equilibrium ratios of the products of the [1,5] hydrogen shift reaction in 2,4,6,7,9-pentamethyl-5-methylene-11,11a-dihydro-12H-naphthacene at 441.05 K. Numbers x,x-d2 denote the deuterium bonding sites.

Conclusions

To conclude, combination of PIMD with higher level methods in the HA proved to be a viable method for accurate calculations of EIEs. The GVE for the derivative of the free energy with respect to the mass allowed us to obtain accurate results at lower temperatures in reasonable time (about 60 times faster than with the TE).
Calculations showed that the effects beyond usual harmonic approximation are important in the computation of the EIE, since they account for up to 30 % of the final value of the reduced free energy of the reactions considered here.
The agreement of theoretically calculated ratios with the data measured in Ref. [2] 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.

Acknowledgements

This research was supported by the Swiss National Science Foundation (Grant No. 200021_124936/1) and by the EPFL.

References

[1] Doering, W.V. and X. Zhao, Effect on kinetics by deuterium in the 1,5-hydrogen shift of a cisoid-locked 1,3(Z)-pentadiene, 2-methyl-10-methylenebicyclo[4.4.0]dec-1-ene: Evidence for tunneling? J. Am. Chem. Soc., 2006. 128(28): p. 9080—9085.
[2] Doering W.V. and E.J. Keliher, Effect of deuterium on the kinetics of 1,5-hydrogen shifts: 5-dideuteriomethylene-2,4,6,7,9-pentamethyl-11,11a-dihydro-12H-naphthace ne. J. Am. Chem. Soc., 2007. 129(9): p. 2488—2495.
[3] Feynman, R.P. and A.R. Hibbs, Quantum Mechanics and Path Integrals. International Series in Pure and Applied Physics, ed. L.I. Schiff. 1965: McGraw-Hill.
[4] Vaníček, J. and W.H. Miller, Efficient estimators for quantum instanton evaluation of the kinetic isotope effects: Application to the intramolecular hydrogen transfer in pentadiene. J. Chem. Phys., 2007. 127(11): p. 114309.
[5] Case, D.A., et al., Amber 10. 2008, University of California: San Francisco.
[6] Frisch, M.J., et al., Gaussian 03, Revision E.01. 2004, Gaussian, Inc.: Wallingford, CT.
[7] Zimmermann T. and Vaníček, J., Path integral evaluation of equilibrium isotope effects. J. Chem. Phys., 2009. 131(2): p. 13.

[1] Isotopologs are molecules which differ by the isotopic substitution of one or more atoms

[2] Isotopomers are molecules which have the same total number of isotopes of its atoms but differ in their positions.



Cherchez ...

- dans tous les Flash informatique
(entre 1986 et 2001: seulement sur les titres et auteurs)
- par mot-clé

Avertissement

Cette page est un article d'une publication de l'EPFL.
Le contenu et certains liens ne sont peut-être plus d'actualité.

Responsabilité

Les articles n'engagent que leurs auteurs, sauf ceux qui concernent de façon évidente des prestations officielles (sous la responsabilité du DIT ou d'autres entités). Toute reproduction, même partielle, n'est autorisée qu'avec l'accord de la rédaction et des auteurs.


Archives sur clé USB

Le Flash informatique ne paraîtra plus. Le dernier numéro est daté de décembre 2013.

Taguage des articles

Depuis 2010, pour aider le lecteur, les articles sont taggués:
  •   tout public
    que vous soyiez utilisateur occasionnel du PC familial, ou bien simplement propriétaire d'un iPhone, lisez l'article marqué tout public, vous y apprendrez plein de choses qui vous permettront de mieux appréhender ces technologies qui envahissent votre quotidien
  •   public averti
    l'article parle de concepts techniques, mais à la portée de toute personne intéressée par les dessous des nouvelles technologies
  •   expert
    le sujet abordé n'intéresse que peu de lecteurs, mais ceux-là seront ravis d'approfondir un thème, d'en savoir plus sur un nouveau langage.