Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2759073

Formats

Article sections

- Abstract
- I. Introduction
- II. Theoretical Background
- III. Computational Details
- IV. Results and Discussions
- V. Conclusions
- Supplementary Material
- References

Authors

Related links

J Am Chem Soc. Author manuscript; available in PMC 2010 October 7.

Published in final edited form as:

PMCID: PMC2759073

NIHMSID: NIHMS146252

The publisher's final edited version of this article is available at J Am Chem Soc

See other articles in PMC that cite the published article.

Primary kinetic isotope effects (KIEs) on a series of carboxylic acid-catalyzed protonation reactions of aryl-substituted α-methoxystyrenes (X-1) to form oxocarbenium ions have been computed using the Kleinert variational second-order perturbation theory (KP2) in the framework of Feynman path integrals (PI) along with the potential energy surface obtained at the B3LYP/6-31+G(d,p) level. Good agreement with the experimental data was obtained, demonstrating that this novel computational approach for computing KIEs of organic reactions is a viable alternative to the traditional method employing Bigeleisen equation and harmonic vibrational frequencies. Although tunneling makes relative small contributions to the lowering of the free energy barriers for the carboxylic acid catalyzed protonation reaction, it is necessary to include tunneling contributions to obtain quantitative estimates of the KIEs. Consideration of anharmonicity can further improve the calculated KIEs for the protonation of substituted α-methoxystyrenes by chloroacetic acid, but for the reactions of the parent and 4-NO_{2} substituted α-methoxystyrene with substituted carboxylic acids, the correction of anharmonicity overestimates the computed KIEs for strong acid catalysts. In agreement with experimental findings, the largest KIEs are found in nearly ergoneutral reactions, Δ*G ^{o}* ≈ 0, where the transition structures are nearly symmetric and the reaction barriers are relatively low. Furthermore, the optimized transition structures are strongly dependent on the free energy for the formation of the carbocation intermediate,

Brønsted acid-catalyzed proton transfer reactions at carbon of conjugated alkenes play an important role in chemical and biological processes.^{1} For example, the bacterial squalene to hopene polycyclization reaction, a process closely related to the biosynthesis of cholesterol, is initiated by the rate-limiting proton transfer from an aspartic acid residue to the 2,3-alkene.^{2}^{-}^{4} An understanding of structure-activity relationship by changing the substituents at the proton donor and acceptor compounds can provide useful insights into the structures of transition state for reactions that generate reactive carbenium ion intermediates both in biological systems and in aqueous solution. In principle, kinetic isotope effects (KIEs) provide a direct probe of the transition state structure;^{5} however, there have been relatively few systematic studies on proton transfer reactions at carbon. The accompanying experimental studies of primary KIEs for the protonation reactions of aryl substituted α-methoxystyrenes **X-1** by a series of substituted acetic acids (**Y**-AcOH)^{6} produce a unique set of experimental data for theoretical analysis to explain the changes in KIE on proton transfer reactions at carbon due to variations in reaction driving force Δ*G ^{o}*. In this article, we use a novel path-integral method based on Kleinert second-order perturbation (KP2) theory to analyze computational results on the same series of reactions to shed light on the origin of the observed trends in isotope effects.

Recently Tsang and Richard reported a fast and general method for determining primary product (PIE) and kinetic isotope effects on proton transfer at carbon in hydroxylic solutions.^{7} That work has been extended to include a broader range of substituents both on the proton donor and on the vinyl ether proton acceptor **X-1**.^{6} It was found that the largest PIEs are observed for the protonation of **4-MeO-1**, which is the most basic and reactive vinyl ether compound in the series. As the basicity decreases in **X-1**, the experimental PIEs show strong dependence on the reaction driving force. These experimental findings are consistent with a Hammond effect on the transition state structure of substituted **X-1** and the strength in acidity of the proton donor.

In this study we perform calculations on the protonation of aryl substituted **X-1** and substituted acetic acids that have been investigated experimentally.^{6}^{,}^{7} Our study is based on Kleinert's second order variational perturbation theory^{8}^{,}^{9} (KP2) for the centroid density of Feynman path integrals^{10} to treat nuclear quantum effects. Our approach goes beyond the traditional method for computing KIE employing the Bigeleisen equation^{11}^{-}^{13} by including contributions from tunneling and vibration anharmonicity. A particularly attractive feature of the KP theory is that the variational perturbation is exponentially and uniformly convergent.^{9} Furthermore, with the instantaneous normal mode approximation,^{14}^{,}^{15} the nuclear quantum effects can be decomposed into anharmonicity and tunneling contributions, providing further insights into the origin of substituent effects on KIE. Although the first order KP theory (KP1), which is equivalent to the variational method independently introduced by Giachett and Tognetti^{16} and by Feynman and Kleinert^{17} (GTFK), has been extensively used in a wide range of applications, higher order KP theory has only been used in limited applications on model systems. We have recently implemented the second and third order Kleinert perturbation theory^{8}^{,}^{9} for molecular systems, in which the path integrals have been integrated analytically.^{14}^{,}^{15}^{,}^{18} This computational procedure is essentially an automated path integral method and integration free since there is no need to carry out discrete Monte Carlo or molecular dynamics simulations. Consequently, it allows the Kleinert perturbation theory (up to the third order) to be conveniently used in practical applications in large molecular systems.

The goals of this study are two folds. First, we aim to provide a theoretical understanding of the trends of the observed primary KIEs on the protonation of **X-1** by different acids.^{6} Secondly, we wish to illustrate the accuracy and applicability of the automated integration-free path-integral (AIF-PI) method^{14} based on KP2 theory for determining kinetic isotope effects of organic reactions using a potential energy surface directly obtained from density functional theory (DFT). We also show that the present KP theory offers an alternative to the traditional approach to compute KIEs for organic reactions using the Bigeleisen and Goeppert-Mayer formalism,^{11}^{-}^{13}^{,}^{19} in which anharmonicity and tunneling contributions are neglected. Of course, the latter effects can be incorporated into the calculation using semiclassical models based on transition state theory, and anharmonic contributions can be analyzed by other separate approaches.^{20}^{,}^{21} Nevertheless, the ability to incorporate these effects in a unified theoretical framework in path integral calculations to efficiently determine KIEs for chemical reactions along with the use of ab initio and DFT calculations can be of great use for organic reactions, making use of Kleinert's variational perturbation theory.^{8}^{,}^{9}

The paper is organized as follows. We first present a brief summary of the theoretical background, focusing on computation of kinetic isotope effects. In section III, the computational details are provided, followed by results and discussions in section IV. In section V, we summarize the major findings from this study.

The most widely used formalism to compute kinetic isotope effects is the Bigeleisen equation, in which the KIE is determined by the use of the harmonic vibrational frequencies, neglecting quantum tunneling:^{11}^{-}^{13}^{,}^{19}

$$\mathrm{KIE}=\left(\frac{{\nu}_{\mathrm{H}}^{\ast}}{{\nu}_{\mathrm{D}}^{\ast}}\right)\left(\prod _{i=1}^{3N-7}\frac{{\nu}_{\mathrm{H}i}^{\ne}}{{\nu}_{\mathrm{D}i}^{\ne}}\prod _{j=1}^{3N-6}\frac{{\nu}_{\mathrm{D}j}}{{\nu}_{\mathrm{H}j}}\right)\left(\prod _{i=1}^{3N-7}\frac{\text{sinh}\left[{u}_{\mathrm{D}i}^{\ne}\right]}{\text{sinh}\left[{u}_{\mathrm{H}i}^{\ne}\right]}\prod _{j=1}^{3N-6}\frac{\text{sinh}\left[{u}_{\mathrm{H}j}\right]}{\text{sinh}\left[{u}_{\mathrm{D}j}\right]}\right)$$

(1)

where *ν*_{Xi} is the *i*th normal mode frequency for an X isotopomer, the superscript ≠ denotes the transition state and the subscript H and D specify the light and heavy isotopes, *u _{i}* =

To go beyond the harmonic approximation and to include quantum tunneling effects in the framework of Feynman path integrals, the Bigeleisen equation can be rewritten as follows:^{18}

$$\mathrm{KIE}=\left(\frac{{\nu}_{\mathrm{H}}^{\ast}}{{\nu}_{\mathrm{D}}^{\ast}}\right)\frac{{e}^{-[{W}_{\mathrm{H}}^{\ne}\left({z}_{\mathrm{H}}^{\ne}\right)-{W}_{\mathrm{D}}^{\ne}\left({z}_{\mathrm{D}}^{\ne}\right)]\u2215{k}_{\mathrm{B}}T}}{{e}^{-[{W}_{\mathrm{H}}^{\mathrm{R}}\left({z}_{\mathrm{H}}^{\mathrm{R}}\right)-{W}_{\mathrm{D}}^{\mathrm{R}}\left({z}_{\mathrm{D}}^{\mathrm{R}}\right)]\u2215{k}_{\mathrm{B}}T}}$$

(2)

where *W*_{X} (*z*_{X}) is the centroid effective potential^{22}^{-}^{24} along the reaction coordinate *z*_{X} calculated from the path integrals for isotopomer X, and ${W}_{\mathrm{X}}^{\ne}\left({z}_{\mathrm{X}}^{\ne}\right)$ and ${W}_{\mathrm{X}}^{\mathrm{R}}\left({z}_{\mathrm{X}}^{\mathrm{R}}\right)$ are its values at the transition state and at the reactant state, respectively. Note that eq 2 reduces to eq 1 when the centroid potential is calculated in the harmonic approximation (and neglecting tunneling effects; Appendix B in Ref. ^{18} provides a proof). Moreover the mass-dependent nature of the centroid effective potential in eq 2 distinguishes it from the classical mechanical potential of mean force in transition state theory, which is mass-independent. Eq 2 is derived based on path-integral quantum transition state theory (PI-QTST),^{22}^{-}^{24} and the reaction coordinate *z*_{X} is defined by the centroid positions of atomic coordinate {* _{k}*}, which is the imaginary time-average ${\stackrel{\u2012}{x}}_{k}=\frac{{k}_{\mathrm{B}}T}{\hslash}{\int}_{0}^{\hslash \u2215{k}_{\mathrm{B}}T}{x}_{k}\left(\tau \right)d\tau $. The mass-dependent property of the centroid potential in the formulation of the path-integral quantum transition-state theory (PI-QTST) has been applied to chemical reactions in condensed phases.

Consequently, the main task in our approach, employing the Kleinert perturbation theory, is to determine the centroid effective potential *W*(*z*) along the reaction coordinate *z*.^{14}^{,}^{15} Although this new approach appears to be very different from the familiar Bigeleisen equation,^{11}^{-}^{13}^{,}^{19} the computational procedure is rather straightforward. As in calculations using eq 1, one first optimizes the reactant state and the transition state structures and determines their harmonic vibrational frequencies. Without including tunneling and anharmonicity effects, eq 1 can be used to estimate KIE at this stage.^{13} In our approach, we further carry out the IRC path calculation to trace the reaction path *z* from the transition sate to the reactant and product state, respectively. In principle, this is required in the use of the Bigeleisen equation for validation of the transition state. The additional calculation needed to construct *W*(*z*) along *z* is to determine the free energy difference between the use of harmonic and anharmonic frequencies plus tunneling contributions by means of centroid path integral method^{10}^{,}^{22}^{-}^{24} using Kleinert perturbation theory.^{8}^{,}^{9}^{,}^{14}^{,}^{15} This is done by using the one-dimensional potential in each normal mode coordinate for the 3*N*-6 normal modes, which is automated in our program.^{14} If the geometry of the transition state is not significantly modified, as in the present case, in the centroid “quantum” reaction coordinate^{22}^{-}^{24} compared with that in the IRC (classical) path, one only needs to perform KP2 calculations at the reactant and transition state geometry, without necessarily following the entire IRC path. The following two paragraphs contain further computational details and the reader who is not interested in these technical issues may skip directly to the next section.

Kleinert's variational perturbation theory provides a systematic approach to determine the path-integral centroid potential *W*(*z*).^{8}^{-}^{10} The derivation for KP theory can be found in Kleinert's textbook^{9} and our computational procedure has been described previously.^{14}^{,}^{15}^{,}^{18} In essence, KP theory systematically builds up anharmonic corrections to a harmonic reference state - not to be confused with the harmonic vibrational frequencies *ν* in eq 1 - defined by a trial angular frequency Ω to yield the *n*th-order perturbation (KPn) approximation to the centroid effective potential $W\left(z\right)\approx {W}_{n}^{\Omega}\left(z\right)$,^{9} which defines the canonical quantum mechanical partition function in path integral representation.^{10} The optimal choice of the trial frequency at a given order of KP expansion and at a particular centroid reaction coordinate is determined by the least-dependence of ${W}_{n}^{\Omega}\left(z\right)$ on Ω, which provides a variational approach to determine its optimal value.^{8}^{,}^{9}^{,}^{14}^{,}^{15}^{,}^{34} The KP theory implemented in the AIF-PI program is an analytical approach allowing the use of high-level *ab initio* or DFT methods to determine the internuclear potential energy function to perform *ab initio* path-integral calculations.^{14}^{,}^{15}^{,}^{18} It would be prohibitively expensive if discretized path-integral Monte-Carlo (PIMC) or centroid molecular dynamics (CMD) simulations^{34} are used for large systems such as the proton transfer reactions (more than 90 degrees of freedom) described here.

An especially attractive feature of KP theory is that if the potential energy surface of the system is expressed as a series of polynomials or Gaussian functions, the path integrals can be explicitly integrated and an analytic expression of the KPn centroid potential ${W}_{n}^{\Omega}\left(z\right)$ can be obtained for a given Ω.^{8}^{,}^{9} However, for multidimensional potentials, the trial frequency becomes a 3*N* × 3*N* matrix for *N* nuclei, which quickly becomes intractable computationally even at the second order perturbation (KP2) level. To render the KP theory computationally feasible for many-body systems, we make instantaneous normal mode (INM) approximation for a given nuclear centroid configuration {* _{k}*}.

$${W}_{n}^{\Omega}\left(z\left\{{\stackrel{\u2012}{x}}_{k}\right\}\right)\approx V\left(\left\{{\stackrel{\u2012}{x}}_{k}\right\}\right)+\sum _{i=1}^{3N}{w}_{n}^{{\Omega}_{i}}\left({q}_{i}\left\{{\stackrel{\u2012}{x}}_{k}\right\}\right)$$

(3)

where *V*({* _{k}*}) is the 3

The molecular structures of the separate reactants (substituted acetic acids and aryl-substituted vinyl ethers), product ion, and the transition-state structures in aqueous solution were optimized using density-functional theory with the inclusion of the polarizable continuum model (PCM)^{31}^{,}^{32} for the treatment of the solvent effects. In all calculations, the hybrid B3LYP functional^{35}^{-}^{37} and the 6-31+G(d,p) basis set were used. Vibrational frequency analyses were performed to confirm the nature of the minimum and saddle points. All the electronic structure calculations were performed using the software package Gaussian 03.^{38} By default, Gaussian 03 uses united atom (UA) model to define the molecular cavity for the PCM model,^{38}^{-}^{40} in which a spherical cavity is placed on each heavy atom encompassing the hydrogen atoms attached to it. To prevent discontinuity of the potential energy surface during transition-state optimizations where the migrating proton is neither completely bonded to the donor (O—H) nor to the acceptor (C—H) atoms, we have chosen to use fixed atomic radii rather than a cavity defined by the isodensity surface. The UAKS radii optimized at the PB0/6-31G(d) level of theory^{39}^{,}^{40} are adopted in the present work except the donor oxygen and the acceptor carbon atoms, whose values have been adjusted to yield the experimental free energy of reaction for *α*-methoxystyrene and chloro-acetic acid. In addition, the hydrogen atoms bonded to the donor and acceptor atoms are explicitly represented in the PCM calculations. All atomic radii used in the present study are given as Supporting Information.

We used eq 2 to determine KIE values in this study. Since the molecular structures of the transition state and the reactant state on the original PCM-B3LYP/6-31+G(d,p) potential energy surface (PES) are similar to those on the centroid PES (their difference is usually in the order of 10^{–2} Å),^{14}^{,}^{15}^{,}^{18} the molecular structures for the transition and reactant states in this paper were not re-optimized on the basis of the centroid potential.

In the present study in path-integral calculations, we express the potential energy surface along each INM coordinate by a twentieth order polynomial that is fitted to results from PCM-B3LYP/6-31+G(d,p) calculations at a step size of 0.1 Å. The centroid potential is computed up to the second order of KP expansion. Thus, the notion for this level of theory is KP2/P20 and P*m* denotes an *m*th order polynomial representation of the original PES of each normal mode coordinate. Through a series of tests,^{14}^{,}^{15}^{,}^{18} we found that KP2/P20 is generally a good choice for the AIF-PI method. The calculated values of the centroid potential are usually within a few percent from the exact values. To compute the anharmonicity contributions to the KIE, we quantized the following five nuclei: the donor oxygen (O1), the protium or deuterium (H1), and the three atoms of the terminal CH_{2} methylene group of **X-1**. The entire system is quantized (more than 30 nuclei) to obtain the imaginary normal mode for computing the tunneling contribution to the KIE. All computations are performed using the AIF-PI program implemented for use with *Mathematica* (Wolfram Research Inc., Mathematica, Versions 5 and 6, Champaign, IL). The analytic expression and program can be found in the Supporting Information of Ref. ^{14} and ^{15}, which is conveniently connected with Gaussian03 calculations, and the program may also be obtained from the authors of this article.

Table 1 lists the bond distances of the transferring proton from the donor (O—H) and the acceptor (C—H) atom, and the experimental and computed PCM free energies of activation Δ*G*^{≠} and free energies of reaction Δ*G ^{o}* (driving force) for the twelve proton transfer reactions, consisting of the entire series of substituted acetic acids and two

Optimized donor (O-H) and acceptor (C-H) bond lengths () at the transition states along with the calculated and experimental free energies of activation (kcal/mol) for the proton transfer reactions between substituted α-methoxystyrenes **...**

Computed total bond orders of the acceptor (O—H) and donor (C—H) bonds *P*_{tot}, the difference in bond order, Δ*P* = *P*(C – H) – *P*(O – H), and in bond lengths (), Δ*R* = *R*(C—H) - *R*(O—H) **...**

The trends of the computed free energies of activation for the proton transfer reactions of the various substituted aryl vinyl ethers (AVE) and substituted acetic acids in water represented by the PCM model are in excellent accord with the experimental data (Table 1), whereas the computed free energies of reaction are somewhat overestimated for the most endoergic reactions.^{6}^{,}^{44} The experimental barrier heights have been obtained from the observed rate constants at 25°C using transition-state theory (TST) and the experimental free energies of reaction are determined from the measured equilibrium constants.^{6}^{,}^{44} For comparison, the root-mean-square-deviation (RMSD) for Δ*G ^{o}* is 3.3 kcal/mol, and the RMSD in the calculated barrier height from the experimental results is 1.2 kcal/mol, in which the largest deviation is about 2 kcal/mol. Overall the computational results are in reasonable agreement with the experimental data, suggesting that the computational method used in the present study is adequate for treating the proton transfer reactions.

The agreement in the trend of free-energy barriers between theory and experiment is mirrored by structural variations and computed bond orders (Table 2). In general, reactions with smaller barriers are accompanied by transition-state structures that are more symmetric with a smaller difference in bond order (Δ*P* = *P*(C – H) – *P*(O – H)) between the acceptor bond and the donor bond (Figure 1). Furthermore, there is a strong linear correlation between the difference in bond order and the asymmetry between the acceptor and donor bond distances (Δ*R* = *R*(C – H) – *R*(O – H)) at the transition state (Figure 2). For the fastest reactions (**4-MeO-1** + ClCH_{2}CO_{2}H and **H-1** + Cl_{2}CHCO_{2}H), Δ*R* is nearly symmetric at 0.003 and 0.002 Å, respectively, and the difference in bond order Δ*P* is 0.157 and 0.156, respectively (Table 2). This is in excellent agreement with the conclusion from analyses of the observed kinetic (product) isotope effects, which also point to more symmetric transition states for reactions with low barriers.^{6} On the other hand, for the slowest reaction of the present study, **4-NO _{2}-1** + CH

Correlation between computed free energy barrier (kcal/mol) and the difference in bond order between the acceptor and donor bond (Δ*P*) at the transition state for the protonation reactions of substituted α-methoxystyrenes by carboxylic **...**

Computed differences in bond length (Δ*R*) and in bond order (Δ*P*) at the transition state for the protonation reactions of substituted α-methoxystyrenes by carboxylic acids.

Computed Hammond shifts in transition structure as illustrated by the difference in acceptor and donor distance from the transferring hydrogen (Δ*R*) with increasing endroergicity of reaction.

The structural variations are reflected by the computed imaginary frequencies at the transition states for the proton transfer reactions. Table 2 shows that the more symmetric and tighter transition states (indicated by small difference in Δ*R*) have the largest values in the imaginary frequency. In fact, a good correlation can be found between the frequency values (either for protium or deuterium) and the differences in bond order of these twelve reactions with an *R*^{2} (correlation coefficient for a linear relationship) of 0.993

The experimental and computed KIE values for the twelve reactions are presented in Table 3. Overall, the agreement between the experiments and calculations is very good with a RMSD of 0.9. The trends observed for structures and bond orders at the transition states as well as the free-energy barriers are clearly reflected in the computed and measured KIEs. Specifically, reactions of the substrate (**X-1**) substituted with a strong electron-donating group (**X** = 4-MeO, and H) or reactions involving a strong carboxylic acids (**Y** = di-Cl, NC-, and Cl), for which the computed transition structures are nearly symmetrical, usually exhibit large KIEs (Table 3). The trend is best depicted by the series of substituted **X-1** with chloroacetic acid (Figure 4). The computed change in KIE is also reasonable for the series of substituted acetic acids (Figure 5) using the parent α-methoxystyrene as the substrate, although there is a leveling effect as the acidity increases in experiments, but not from computaton.^{6} One possible reason for this small difference is the use of a single fractionation factor for different substituted acids derived from acetic acid to convert the experimental PIEs into KIEs.^{7}^{,}^{45} However, the computational results include the difference in fractionation factors among different acids (see below). For example, if a smaller fractionation factor is used for the weaker acids in the experimental data, the agreement with theory can be significantly improved.

Variations of the computed KIEs using the Bigeleishen equation with harmonic frequencies (green) along with the inclusion of anharmonicity (blue) and tunneling (red) for the reactions of α-methoxystyrene with various substituted acetic acids ( **...**

Computed and experimental primary kinetic isotope effects (KIE), along with the KIE determined using the Bigeleisen equation with harmonic vibrational frequencies, the anharmonicity correction factor and tunneling contributions, for the protium/deuterium **...**

Figure 6 shows the temperature-dependence of the KIE for the reactions of *α*-methoxystyrene (**H-1**) and methoxyacetic acid computed using KP2/P20. A slope of $-({E}_{a}^{\mathrm{H}}-{E}_{a}^{\mathrm{D}})\u2215R=518\phantom{\rule{thickmathspace}{0ex}}\mathrm{K}$, where *R* is the gas constant, and an intercept of ln(*A*^{H} / *A*^{D}) = 0.128 are obtained from Figure 6. Both values are in reasonable agreement with the experimental results, which are 595 ± 30 K and 0.02 ± 0.10, respectively. This yields a computed difference of ${E}_{a}^{\mathrm{H}}-{E}_{a}^{\mathrm{D}}=-1.03$ kcal/mol for the protonation of **H-1** by H and D, which may be compared with the experimental value of −1.18 kcal/mol. The slightly underestimated ${E}_{a}^{\mathrm{H}}-{E}_{a}^{\mathrm{D}}$ difference is reflected by a smaller calculated KIE (6.45) than that from experiments (6.8). The Y-intercept in Figure 6 gives a value of 1.14 for the ratio *A*^{H} / *A*^{D}, but the experimental value is closer to unity (1.02). The computed temperature dependence in the Arrehnius-type plot indicates that there is a small amount of hydrogen tunneling, but the extrapolation of the experimental temperature dependence of KIE indicates that there is little tunneling.^{6}^{,}^{7}^{,}^{44} Nevertheless, we note that the computed tunneling effects contributes no more than 0.1 kcal/mol to the overall thermal activation barrier of 18.4 kcal/mol (17.5 kcal/mol from experiment) for the protonation of **H-1** by methoxyacetic acid, and these proton transfer reactions are dominated by over barrier processes. Figures 5 and and66 (see also Table 3) show that the computed KIEs are dominated by contributions from harmonic vibrations, whereas anharmonicity and tunneling contributions are relatively small, but they make important contributions to the calculation of the absolute KIE values.^{5}^{,}^{21}

In the present integration-free path-integral method,^{14}^{,}^{15} we make use of independent normal mode coordinates in the calculations of the centroid potential of path integrals. Although this assumption neglects correlation between different normal mode coordinates, the error introduced appears to be small from previous studies^{14}^{,}^{15}^{,}^{18} and for the present proton transfer reactions as evidenced by the good agreement in the computed free-energy barriers and KIE values with experiments.^{6} Importantly, it allows us to separate the overall quantum effects into contributions from harmonic vibrations,^{21} anharmonicity, and quantum tunneling to provide further insight into the origin of the dependence of KIE on substituent effects, both on the substrate and the proton donor. In this analysis, the total KIE in eq 2 can be decomposed as follows:

$$\mathrm{KIE}={\eta}_{har}^{KIE}{h}_{anh}^{KIE}{\kappa}_{tun}$$

(4)

where ${\eta}_{har}^{KIE}$ is the KIE determined using the Bigeleisen equation,^{11}^{-}^{13} eq 1, with harmonic vibrational frequencies, ${h}_{anh}^{KIE}$ is the anharmonicity correction factor for the KIE from all modes with real (positive) frequencies as determined by the KP2 theory of the centroid path integrals using the exact potentials from PCM-B3LYP/6-31+G(d,p) calculations, and *κ _{tun}* represents tunneling contributions to the computed KIE value, corresponding to the imaginary normal modes. Refs.

The computed KIEs using the Bigeleisen equation (eq 1) for the three series of protonation reactions of **X-1** and **Y**-AcOH(D) are shown in Table 3 and Figures 4 and and5.5. The overall trends of the computed “harmonic” KIEs are in accord with the corresponding experimental data;^{6} however, they are underestimated by more than 30% for the fast reactions and by 10-20% for the slower processes. Quantum mechanical tunneling is undoubtedly important in computing kinetic isotope effects,^{5}^{,}^{20}^{,}^{21} and inclusion of tunneling contributions by multiplying the tunneling factor in Table 3 significantly improves the computed KIE for the slower protonation reactions of substituted α-methoxystyrenes, reducing the errors to about 5% or less, but errors remain high for the faster reactions at nearly 20% of the experimental results. These reactions include those using very strong acids as the proton donor and electron-rich acceptors (**4-MeO-1**). In these cases, it is important to go beyond the Bigeleisen approach to obtain more accurate KIEs.

Figure 4 shows that the KIE determined using harmonic frequencies (green curve) decreases with increasing electron-withdrawing power of the substituent on the aryl ring. This may be attributed to the effects of charge polarization of the vinyl double bond with a reduced bond order. Consequently, there is somewhat reduced zero-point effect from the substrate **X-1** in going from the reactant state to the transition state. This is consistent with the trend in Figure 5, in which the KIEs computed using harmonic frequencies are similar. Figure 5 further shows that there is little variation from harmonic frequencies with different donor acidities and this suggests that the amount of zero-point energy from the donor hydroxyl bond that is maintained in the symmetric stretch vibration is minimal.^{46} However, the present results do not address this issue directly, because there are other factors, such as the change in the bond order of the substrate X-1 and donor Y-AcOH.

The quantum tunneling effect has a strong dependence on the symmetry of the donor and acceptor bond distance from the transferring proton (deuteron) at the transition state (by definition, tunneling only occurs near the transition state). Examination of Tables 2 and and33 shows that the ratio of tunneling for the H/D transfer reactions, *i.e.*, the tunneling factor in eq 4 (*κ _{tun}*), has an inverse exponential correlation with the difference in bond lengths, Δ

Tables 3 and and44 list the anharmonicity correction factors for the three series of reactions between substituted **X-1** and substituted acetic acid (**Y**-AcOH) in water. The percent contribution to the total computed KIE (*k ^{H}* /

$${h}_{anh}^{KIE}=\frac{{h}_{anh}^{TS}}{{h}_{anh}^{RS}}$$

(5)

where ${h}_{anh}^{RS}$ and ${h}_{anh}^{TS}$ specify the isotope effects on anharmonicity in the reactant state (RS) and the transition state (TS), respectively. ${h}_{anh}^{RS}$ accounts for the isotope effects on anharmonicity of the acidity constant for exchange of H and D at the carboxylic acid, *i.e.*, the anharmonic contributions to the factor 1/ Φ_{AL}, which is the equilibrium constant for the exchange of H and D by the carboxylic acid in 50/50 mixture of H_{2}O/D_{2}O.^{6} It is the ratio of the Boltzmann factors for the H and D vibrational energies computed using harmonic frequencies and with the real potential energy surface that includes anharmonicity in a given acid. Thus, ${h}_{anh}^{RS}$ can be calculated by

$${h}_{anh}^{RS}=\frac{{e}^{-\Delta {W}_{har\to anh}^{RS}\left(H\right)\u2215{k}_{B}T}}{{e}^{-\Delta {W}_{har\to anh}^{RS}\left(D\right)\u2215{k}_{B}T}}$$

(6)

where $\Delta {W}_{har\to anh}^{RS}$ is the change in quantum energy (mostly zero-point energy and thermal contributions) due to anharmonicity at the carboxylic acid. In other words, the experimental fractionation factor defined by eq 11 in Ref. ^{6} is related to that computed using harmonic frequencies multiplied by ${h}_{anh}^{RS}$, *i.e.*, $(1\u2215{\Phi}_{\mathrm{AL}}^{\text{exp}})={h}_{anh}^{\mathrm{AL}}(1\u2215{\Phi}_{\mathrm{AL}}^{\text{har}})$, where AL *RS*.^{6} The exchange for H and D at the TS, ${h}_{anh}^{TS}$, may be similarly defined and can be interpreted as the anharmonic contribution to the residual zero-point energy maintained at the transition state for the transferring H and D.^{6}

Computed isotope effects on anharmonicity in the reactant state and the transition state along with the total anhamonic correction factors the proton transfer reactions from substituted acetic acids to aryl-substituted α-methoxystyrenes.

Table 4 lists the total anharmonicity correction factors (Table 3) and the isotope effects on anharmonicity at the Brønsted acid and the transition state. Examining the fourth column of Table 4, we immediately notice that the isotope effects on anharmonicity are close to unity and nearly constant for the exchange of H and D at the transition state for all reactions considered. There is no clear dependence either on the substituents of **X-1** or on the acidity strength of the proton donor in the small variations of a few percent (Table 3). However, there are noticeable isotope effects on anhamonicity for the exchange of H and D at the substituted carboxylic acids, with ${h}_{anh}^{RS}$ deviating from unity by 10-22% (except the parent reaction of **H-1** with acetic acid). Furthermore, one can identify that the change in ${h}_{anh}^{RS}$ is closely correlated with the change in the electron-withdrawing power of the acid substituent (*i.e.*, its acidity). This suggests that the anharmonic effect on the computed KIEs is mainly due to its modification of the H/D enrichment (Φ_{AH}) of the Brønsted acid by using eq 1 with harmonic vibrational frequencies.

We have identified three key normal modes that have the largest anharmonic contributions (normal mode #3, #5, and #6), and they are all related to the protium/deuterium (H/D) being transferred (Figure 8). In Figure 9, we illustrate the difference in quantum mechanical effects, mostly zero-point-energy (ZPE) effects, between **Y**-AcO(H) and **Y**-AcO(D). The differences in the total quantum harmonic vibrational energy or harmonic centroid potential energy are shown in green, decomposed into the three specific modes. The corresponding anharmonic contributions are illustrated in red. The sum of harmonic and anharmonic corrections is the total quantum energy difference from the centroid potential between the two isotopic acids. The trend in Figure 9 follows the increase in donor acidity, resulting in a decrease in bond strength and force constants involving the H/D modes. This is reflected by the decreasing trend of ZPE from harmonic modes (green bars in Figure 9). However, as the force constants decrease, (the potential energy surface becomes flatter) and the deviation from the harmonic approximation becomes greater, leading to larger anharmonic contributions. This is consistent with the increasingly large anharmonic energy corrections in Figure 9 (red). This is reflected by the observation and trends of anharmonic contributions to the computed KIE in Table 4.

Schematic illustration of the three normal modes identified to have the largest anhamonic effects as electron-withdrawing substituents are introduced. All three modes are associated with the protium/deuterium as depicted by chloroacetic acid.

Finally, in Figure 10, we illustrate the variation in KIE as a function of the thermodynamic driving force, Δ*G ^{o}*. The experimental data cover a large range of 16 kcal/mol, which are compared with the computational data for the three series of different substitutions on

Primary kinetic isotope effects (KIEs) on a series of carboxylic acid-catalyzed protonation reactions of aryl-substituted α-methoxystyrenes (**X-1**) to form oxocarbenium ions have been computed using the Kleinert variational second-order perturbation theory (KP2) in the framework of Feynman path integrals (PI) and the B3LYP/6-31+G(d,p) potential surface coupled with the polarizable continuum model to represent aqueous solution. Reasonable agreement with the experimental data was obtained, demonstrating that this novel computational approach for computing KIEs of organic reactions is a viable alternative to go beyond the traditional method employing Bigeleisen equation, where anharmonicity and tunneling contributions are important. We found that the dominant factor contributing to the observed KIE for the protonation reactions of substituted α-methoxystyrenes is due to the change in zero-point energy from the reactant to the transition state, and the trend of the computed KIEs using Bigeleisen equation is in accord with experiment. However, it is necessary to include tunneling contributions to obtain quantitative estimates of the KIEs. Without tunneling, the RMSD error in the computed KIEs using Bigeleisen equation is 23% in comparison with the experimental data, and this is reduced to 15% with tunneling contributions. Consideration of anharmonicity can further improve the calculated KIEs; for example, for the protonation of substituted **X-1** by a common acid, chloroacetic acid, the anharmonicity corrected KIEs are in quantitative agreement with experiment (Figure 4). On the other hand, for the reactions of α-methoxystyrene and **4-NO _{2}-1** with substituted carboxylic acid, the correction of anharmonicity leads to overestimation of the computed KIEs for strong acid catalysts, but there is a noticeable leveling effect in the observed KIEs as the acidity increases (Figure 5).

The results for experiments and calculations are each consistent with the conclusion that proton transfer from Brønsted acids to **X-1** proceeds mainly by an activation-limited reaction, where the proton passes over an energy barrier. However, extrapolation of a linear Arrhenius type relationship of the product isotope effects (PIEs) for protonation of **H-1** by methoxyacetic acid to 1/T = 0 gives a limiting PIE of 1.0 at infinite temperature. This result provides no evidence for the existence of a tunneling pathway at room temperature, while the tunneling correction results in an increase in the isotope effect over that calculated for proton transfer over a reaction barrier (*i.e.*, without tunneling), ranging from as little as 1% to as large as 26%. We conclude that either these Arrhenius type relationships are unable to detect relatively small amounts of tunneling for these room temperature proton transfer reactions, or that these calculations overestimate the required tunneling correction. In agreement with the experimental findings in the accompanying paper, the largest KIEs are found in reactions with Δ*G ^{o}* ≈ 0, where the transition structures are nearly symmetric and the reaction barriers are relatively low. Furthermore, we found that the optimized transition structures are strongly dependent on the free energy for the formation of the carbocation intermediate, Δ

This work was partially supported by the National Institutes of Health under grants GM46736 (J.G.) and GM39754 (J.P.R.).

1. Wendt KU, Schulz GE, Corey EJ, Liu DR. Angew. Chem. Int. Ed. 2000;39:2812. [PubMed]

2. Wendt KU, Poralla K, Schulz GE. Science. 1997;277:1811. [PubMed]

3. Cheng J, Hoshino T. Org. Biomol. Chem. 2009;7:1689. [PubMed]

4. Rajamani R, Gao J. J. Am. Chem. Soc. 2003;125:12768. [PubMed]

5. Kohen A, Limbach H-H. Isotope effects in chemistry and biology. Taylor & Francis; Boca Raton: 2006. p. 1074.

6. Tsang W-Y, Richard JP. accompanying paper.

7. Tsang W-Y, Richard JP. J. Am. Chem. Soc. 2007;129:10330. [PMC free article] [PubMed]

8. Kleinert H. Phys. Lett. A. 1993;173:332.

9. Kleinert H. Path integrals in quantum mechanics, statistics, polymer physics, and financial markets. 3rd ed. World Scientific; Singapore; River Edge, NJ: 2004. p. xxvi.p. 1468. For the quantum mechanical integral equation, see Section 1.9; For the variational perturbation theory, see Chapters 3 and 5.

10. Feynman RP, Hibbs AR. Quantum mechanics and path integrals. McGraw-Hill; New York, NY: 1965. p. 365.

11. Bigeleisen J. J. Chem. Phys. 1949;17:675.

12. Bigeleisen J. J. Chem. Phys. 1955;23:2264.

13. Schaad LJ, Bytautas L, Houk KN. Can. J. Chem. 1999;77:875.

14. Wong K-Y, Gao J. J. Chem. Phys. 2007;127:211103. [PMC free article] [PubMed]

15. Wong K-Y, Gao J. J. Chem. Theory Comput. 2008;4:1409. [PMC free article] [PubMed]

16. Giachetti R, Tognetti V. Phys. Rev. Lett. 1985;55:912. [PubMed]

17. Feynman RP, Kleinert H. Phys. Rev. A. 1986;34:5080. [PubMed]

18. Wong K-Y. Ph.D. Thesis. University of Minnesota (USA); 2008.

19. Wolfsberg M. In: Isotope Effects in Chemistry and Biology. Kohen A, Limbach H-H, editors. Taylor & Francis; Boca Raton: 2006. pp. 89–117.

20. Fernandez-Ramos A, Miller JA, Klippenstein SJ, Truhlar DG. Chem. Rev. 2006;106:4518. [PubMed]

21. Pu J, Gao J, Truhlar DG. Chem. Rev. 2006;106:3140. [PMC free article] [PubMed]

22. Gillan MJ. Phys. Rev. Lett. 1987;58:563. [PubMed]

23. Cao J, Voth GA. J. Chem. Phys. 1994;101:6168.

24. Voth GA. Adv. Chem. Phys. 1996;93:135.

25. Voth GA, Chandler D, Miller WH. J. Chem. Phys. 1989;91:7749.

26. Jang S, Schwieters CD, Voth GA. J. Phys. Chem. A. 1999;103:9527.

27. Jang S, Voth GA. J. Chem. Phys. 2000;112:8747–8757. Erratum: 2001, 114, 1944.

28. Gao J, Wong K-Y, Major DT. J. Comput. Chem. 2008;29:514. [PMC free article] [PubMed]

29. Wang M, Lu Z, Yang W. J. Chem. Phys. 2006;124:124516. [PubMed]

30. Wang Q, Hammes-Schiffer S. J. Chem. Phys. 2006;125:184102. [PubMed]

31. Miertus S, Scrocco E, Tomasi J. Chem. Phys. 1981;55:117.

32. Cammi R, Tomasi J. J. Comput. Chem. 1995;16:1449.

33. Fukui K. Acc. Chem. Res. 1981;14:363.

34. Poulsen JA, Nyman G, Rossky PJ. J. Chem. Theory Comput. 2006;2:1482. [PubMed]

35. Lee C, Yang W, Parr RG. Phys. Rev. B. 1988;37:785. [PubMed]

36. Becke AD. J. Chem. Phys. 1993;98:1372.

37. Stephens PJ, Devlin FJ, Chabalowski CF, Frisch MJ. J. Phys. Chem. 1994;98:11623.

38. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Montgomery JA, Jr., Vreven T, Kudin KN, Burant JC, Millam JM, Iyengar SS, Tomasi J, Barone V, Mennucci B, Cossi M, Scalmani G, Rega N, Petersson GA, Nakatsuji H, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Klene M, Li X, Knox JE, Hratchian HP, Cross JB, Bakken V, Adamo C, Jaramillo J, Gomperts R, Stratmann RE, Yazyev O, Austin AJ, Cammi R, Pomelli C, Ochterski JW, Ayala PY, Morokuma K, Voth GA, Salvador P, Dannenberg JJ, Zakrzewski VG, Dapprich S, Daniels AD, Strain MC, Farkas O, Malick DK, Rabuck AD, Raghavachari K, Foresman JB, Ortiz JV, Cui Q, Baboul AG, Clifford S, Cioslowski J, Stefanov BB, Liu G, Liashenko A, Piskorz P, Komaromi I, Martin RL, Fox DJ, Keith T, Al-Laham MA, Peng CY, Nanayakkara A, Challacombe M, Gill PMW, Johnson B, Chen W, Wong MW, Gonzalez C, Pople JA. Gaussian 03, Revision C.02. Gaussian, Inc.; Wallingford CT: 2004.

39. Barone V, Cossi M, Tomasi J. J. Chem. Phys. 1997;107:3210.

40. Takano Y, Houk KN. J. Chem. Theory Comput. 2005;1:71. [PubMed]

41. Wiberg KB. Tetrahedron. 1968;24:1083–96.

42. Glendening ED, Weinhold F. J. Comput. Chem. 1998;19:610–627.

43. Glendening ED, Badenhoop JK, Weinhold F. J. Comput. Chem. 1998;19:628.

44. Richard JP, Williams KB. J. Am. Chem. Soc. 2007;129:6952. [PMC free article] [PubMed]

45. Jarret RM, Saunders M. J. Am. Chem. Soc. 1985;107:2648.

46. Westheimer FH. Chem. Rev. 1961;61:265–273.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |