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

**|**HHS Author Manuscripts**|**PMC2992356

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Theoretical Background
- 3. Simulation Details
- 4. Results & Discussion
- 5. Concluding Remarks
- Supplementary Material
- References

Authors

Related links

J Chem Theory Comput. Author manuscript; available in PMC 2011 September 14.

Published in final edited form as:

J Chem Theory Comput. 2010 September 14; 6(9): 2566–2580.

doi: 10.1021/ct900579kPMCID: PMC2992356

NIHMSID: NIHMS227183

Kim F. Wong,^{a,}^{†} Jason L. Sonnenberg,^{b,}^{‡} Francesco Paesani,^{a,}^{±} Takeshi Yamamoto,^{c} Jiří Vaníček,^{d,}^{§} Wei Zhang,^{e} H. Bernhard Schlegel,^{b,}^{*} David A. Case,^{e} Thomas E. Cheatham, III,^{f} William H. Miller,^{d,}^{*} and Gregory A. Voth^{a,}^{*}

See other articles in PMC that cite the published article.

The rates of intramolecular proton transfer are calculated on a full-dimensional reactive electronic potential energy surface that incorporates high level ab initio calculations along the reaction path and by using classical Transition State theory, Path-Integral Quantum Transition State Theory, and the Quantum Instanton approach. The specific example problem studied is malonaldehyde. Estimates of the kinetic isotope effect using the latter two methods are found to be in reasonable agreement with each other. Improvements and extensions of this practical, yet chemically accurate framework for the calculations of quantized, reactive dynamics are also discussed.

The breaking and formation of chemical bonds is fundamental to chemistry. While molecular dynamics (MD) simulations of large-scale assemblies for hundreds of nanoseconds and even for a few microseconds are permissible using current terascale and emergent petascale high-performance computing infrastructures,^{1}^{-}^{3} the majority of these scientific applications involving conformational sampling or molecular association processes cannot model chemical reactions. Much of current force field research focuses on the development of accurate and reliable interactions for describing protein structures, solvation energies and hybrid bio-inorganic material properties.^{4}^{-}^{10} Designing a chemically reactive potential energy surface (PES), nevertheless, is not a new idea.^{11} During the early half of the twentieth century, chemical physics pioneers, such as Eyring,^{12} Evans and Polanyi,^{13}^{-}^{15} proposed procedures for forming surfaces to describe diverse systems involving ionic, S_{N}2 and Diels-Alder reactions. The common element among these methodologies is the description of the system within a matrix representation for the Hamiltonian

(1)

(2)

Here the system evolves on the lowest energy eigenstate, *ε*_{0}, of the Hamiltonian as a superposition state of multiple nuclear-electronic configurations.^{16} Variations of this technique differ primarily in the treatments of the *H _{ii}* and

Although the empirical approach is a computationally practical strategy for modeling chemistry, a general framework for constructing reactive potential energy surfaces based on *first principles*, ab initio, information is desirable. Empirical procedures typically involve a high startup cost in the form of parameter fitting and calibrations that in some instances may be system specific. On the other hand, even with current advances in computation infrastructure (both from an algorithms and hardware perspective), direct ab initio MD approaches are limited from small to moderate size systems for a few hundreds of ps of conformational sampling.^{19}^{-}^{24} These timescales cannot capture the physics of many biologically relevant reactions, typically occurring in the microseconds or longer timescales. Furthermore, the computation time required for adequately sampling phase space to obtain converged thermodynamic observables precludes the use of direct ab initio MD approaches for complex systems consisting of 10,000 atoms or more. A general hybrid simulation framework combining the computational advantages of MD simulation with the accuracy of electronic structure methods is appealing.

In the present paper, we describe an implementation of such a general framework, the distributed Gaussian (DG) approach,^{25}^{-}^{27} for developing accurate ab initio-based potential energy surfaces capable of modeling chemical reactions. The purpose of this paper is to describe the integration, within the Amber^{28} biosimulations suite, of the DG-EVB surface generation methodology with molecular dynamics, Path Integrals,^{29}^{-}^{31} and Quantum Instanton^{32}^{-}^{35} approaches for the computation of thermal rate constants. This methodology development is tested on a well established system, the prototype intramolecular proton transfer reaction in malonaldehyde coupled to a thermal bath. Despite the presumed simplicity of malonaldehyde, it is only recently that full-dimensional quantum calculations of tunneling splitting based on accurate ab initio surfaces were possible.^{36}^{-}^{40} Previous works utilized reduced-dimensional approximations or less accurate potentials that may not fully capture the physics of the reaction.^{41}^{-}^{47} Malonaldehyde, thus, provides a well studied yet challenging test for assessing the diverse components of the DG-EVB methodology. Integrating both molecular mechanical and ab initio elements, DG-EVB is similar in strategy to Truhlar's recently developed multi-configuration molecular mechanics (MCMM) method.^{48}^{-}^{52}

The remainder of this paper is organized as follows: To show the relationship between the various methods, the theoretical backgrounds for the DG method, the path integral approach and the QI method are briefly reviewed. The simulation details are described in Section 3 followed by the Results and Discussion in Section 4. Lastly, we conclude with an assessment of the method and suggestions for future improvements.

This section briefly describes the theoretical basis to our algorithmic development for constructing a reactive PES from ab initio information and the incorporation of nuclear quantum effects in molecular dynamics. Dynamical methods for computing thermal rate constants and kinetic isotope effects are described. Our test system of a single malonaldehyde molecule coupled to a thermal reservoir is modeled within the canonical ensemble. Within the canonical ensemble, thermal rate constants are well defined quantities.^{53} Under favorable circumstances, tunneling splittings can also be extracted from the imaginary time correlation function from PIMD simulations, but in the present case reliable error bounds on the tunneling splitting could not be obtained using the maximum entropy approach.

A number of methods for building an ab initio-derived PES using Eq. (1) has been proposed recently. Among these are Chang-Miller,^{54}^{-}^{55} Minichino-Voth,^{56} multi-configuration molecular mechanics (MCMM)^{48}^{-}^{52} and our distributed Gaussian EVB (DG-EVB).^{25}^{-}^{27} Since the distributed Gaussian approach is related to the Chang-Miller prescription, we concisely describe Chang-Miller below to motivate the subsequent method development. The Chang-Miller approach attempts to construct an accurate reactive potential energy surface by fitting a superposition of reactant and product configurations using a generalized Gaussian form for the coupling term

(3)

The diagonal elements describing the reactant state (RS) and product state (PS) are approximated using conventional, non-reactive force fields (FF). This choice was motivated by the expectation that the development of classical force fields will continually improve towards a state where they are accurate *enough*. Formulated to reproduce the ab initio energy *ε*_{Ψ}, gradient and Hessian at the transition state geometry, **q**_{TS-}, the parameters *A* (a scalar), **B** (a vector) and (a matrix) take the analytical form

(4)

(5)

(6)

When the system configuration deviates far from the transition state structure (i.e., for large Δ**q**) and the matrix contains one or more negative frequencies,
diverges. Although refinements are available for controlling the asymptotic behavior of the Chang-Miller approach,^{55} simply recasting Eq. (3) in terms of a quadratic polynomial times a spherical Gaussian

(7)

keeps the coupling element bounded at the asymptotes.^{25} The scalar *A* and vector **B** parameters are identical to those in the Chang-Miller approach, while the matrix assumes a modified form

(8)

Note that **Ĩ** is the identity matrix and *α* is a parameter related to the Gaussian width. The distributed Gaussian approach generalizes the above polynomial times a Gaussian prescription [Eq. (7)] to utilize ab initio information not only at the transition state geometry but also at other points on the potential energy surface. Here,
(**q**) is approximated as an expansion in a set of *distributed Gaussians* centered on a set of molecular configurations **q**_{K}

(9)

(10)

(11)

(12)

Δq_{K} = q − q_{K}

(13)

where *N*_{cfg} is the number of ab initio data points used for the fitting, *N*_{dim} is the number of system coordinates, *g*(**q**, **q*** _{K}*,

(14)

that can be solved using singular value decomposition or by an iterative procedure, such as GMRES (generalized minimal residual method).^{57}^{-}^{60} When derivatives for the coupling terms are available, this information can also be utilized for the fitting, i.e.,

(15)

(16)

The **F** column vector stores the terms
(**q*** _{L}*),
and
evaluated at the

For a symmetric 2×2 Hamiltonian matrix, the analytical expression for the coupling and derivatives are given by

(17)

(18)

(19)

where *ε*_{0} is the lowest eigenvalue of the matrix and *H*_{11} and *H*_{22} are the reactant and product valence bond states. *H*_{11} and *H*_{22} can be described by a force field potential (as in the spirit of Chang-Miller^{54}) or as a Taylor series expansion about the respective ab initio minimum.^{25} The key idea is to determine
(**q**) such that the resulting *ε*_{0} surface approximates the ab initio surface, i.e., *ε*_{0} = *ε*_{Ψ}. For the above two-state system, the coupling and corresponding derivatives can be evaluated directly from Eqs. (17)-(19) using the ab initio energies (*ε*_{Ψ}), gradients (*ε*_{Ψ}/**q**) and Hessian (^{2}ε_{Ψ}*/***q**^{2}) data in conjunction with the corresponding *H _{ii}* values and derivatives. This procedure provides coupling values and derivatives for the

When the reactive process involves particles whose thermal de Broglie wavelength is comparable to the characteristic scale of the potential, nuclear quantum effects such as tunneling and zero-point motion are important for describing the chemical rate.^{61}^{-}^{63} Feynman's Path Integral (PI)^{29}^{-}^{31} formalism of quantum mechanics provides a general framework for describing equilibrium (i.e., thermodynamic and structural) properties of a quantum many-body system. Briefly, the quantum partition function for the canonical (NVT) ensemble can be written as

(20)

where

(21)

In this form, *Q* is isomorphic to the classical configurational partition function of *N* ring polymers, each comprised of *P* particles.^{64} Each particle within the polymer interacts harmonically with neighbors along the cyclic path (first term) and with other particles within the same *s*-indexed PI slice via *V*(**q**^{(s)}) / *P* (second term). By introducing a set of Gaussian integrals into Eq. (20), one sees that the quantum partition function

(22)

can be evaluated from molecular dynamics based on the Newtonian equations of motion derived from a fictitious classical Hamiltonian of the form

(23)

This is possible because the inserted fictitious momenta are uncoupled and thus can be integrated analytically such that the prefactor Λ can be defined to give back the correct quantum partition function in Eq. (20). The purpose of the fictitious Hamiltonian is simply to provide a computationally practical scheme for evaluating the total phase space integral. Finally, to recover *Q* in the NVT ensemble, we also need to couple the system to a thermostat so as to ensure that the sampled distribution is indeed canonical.^{65}

While the path integral formulation is general for any inter-nuclear potential, the results will depend on the accuracy of the electronic surface. In our study, the DG-EVB potential enters into PIMD through the second term of Eq. (21). The PI results presented here are computed from the normal mode representation of PIMD,^{66} where the thermostat is a set of Nosé-Hoover chains.^{67}

The computable quantities connecting theory and simulation to experimental measurements are traditionally referred to as *observables*. For a reactive process, the *naturally* measurable quantity is the chemical rate. Transition state theory (TST)^{53}^{,}^{68}^{-}^{70} provides a simple framework for calculating the rate from molecular dynamics trajectories. In TST, the rate along a *reaction path ξ* (**q**) [synonymous with the *reaction coordinate* (RC)] is given by

(24)

where *dξ/dt* is the *dynamical frequency factor* and *ρ* is the normalized density of sampled RC values

(25)

We note that the *intrinsic reaction coordinate* (minimum-energy reaction path) is one well-defined prescription for *ξ* (**q**), although any RC (subject to the limitations of TST) may be employed. In the above equation, *β* = 1/*k _{B}T*,

(26)

where _{ξ≠} denotes the conditional average computed at the dividing surface *ξ*^{≠}. Both the average distribution of RC values in Eq. (25) and the gradient of the RC in Eq. (26) are readily computable from MD using umbrella sampling techniques, for example.

For pedagogical purposes, it is convenient to recast the average distribution function

(27)

in the perspective of the *potential of mean force* (PMF)^{71}

(28)

where *ξ*^{0} and *w***(***ξ*^{0}**)** are constants that are typically chosen to reflect initial conditions of the chemical system. For example, note that *ρ***(***ξ***)** and **(***ξ***)** in the above equations differ only by a normalization factor. This factor, within transition state theory, is chosen to normalize the initial distribution in the reactant state region.

The PMF relates the probability of sampling along the RC to a *free energy* profile. Higher probabilities are associated with relatively low free energy values compared to regions of lower probabilities. At the bottleneck (i.e., at min [**(***ξ***)**]), the PMF corresponds to the barrier for a process under observation. The PMF, therefore, provides an intuitive *free energy surface* perspective to chemical dynamics. In terms of *w***(***ξ***)**, the TST normalized density factor can be rewritten as

(29)

where the limits of integration in the denominator is over the reactant state region.

Another experimental measure of reactive chemical dynamics, capable of providing insights into nuclear quantum effects, is the kinetic isotope effect (KIE). The KIE is defined as the ratio of the reaction rate, *k*, involving the system with a lighter isotope (*l*) compared to the rate involving the system with a heavier (*h*) isotopic substitution

(30)

When the isotopic substitution involves a chemical bond pertaining to the rate limiting step, the KIE is characterized as *primary*. If the substitution does not directly involve chemical bonds that are broken and formed, the KIE is referred as *secondary*. This paper only addresses the primary KIE for the intramolecular proton transfer reaction in malonaldehyde; however, the methodology presented here can be applied to estimating the secondary KIE as well.

As described in the previous section, the inclusion of zero-point motion and nuclear tunneling are important for calculating the rate of proton transfer. These nuclear quantum effects can be taken into account by using a *quantum* analog of TST based on the path integral formalism^{72}

(31)

where *dξ _{c}/dt* is the dynamical frequency factor and

(32)

where *β* = 1/*k _{B}T*,

(33)

Furthermore, the centroid density can be recast in terms of a *quantum* PMF *w***(***ξ _{c}*

(34)

where the integration limits are again over the RS region.

From Eq. (31), it is possible to compute the KIE for proton transfer in malonaldehyde by directly forming the ratio of the rates of reaction, i.e.,
, where (H) denotes the system with the hydrogen isotope and (D) indicates the deuterium isotopic substitution. Alternatively, one can estimate the ratio of the hydrogen and deuterium centroid densities using thermodynamic integration (TI) with respect to mass^{73}^{-}^{74}

(35)

where the parameter *λ* interpolates between the masses of the hydrogen and deuterium isotopes as

(36)

Substituting the centroid density expression [Eq. (32)] into the logarithmic derivative, we obtain

(37)

where _{RS} denotes the canonical average over the reactant state region,
denotes the conditional average with the system constrained to the dividing surface
and *d*Φ/*dλ* can be computed from the PI virial-like estimator^{65}^{,}^{75}

(38)

All elements in the above equation are previously defined. While two separate trajectories are required for each value of *λ*, the savings in computation (compared to explicit construction of the PMFs) can be substantial if one needs to performed multiple umbrella sampling trajectories in order to map out the entire RC range due to an intrinsic high free energy barrier. *d*Φ/*dλ* is typically a smooth, slowly-varying function and thus can be numerically integrated using the simple Trapezoidal Rule along uniformly spaced *λ* points. Because the isotopic ratio of the dynamical frequency factors only involve conditional averages with the system constrained to
, the computation is not expensive and thus it is not necessary to consider TI for this factor.

An alternative to PI-QTST for incorporating nuclear quantum effects in thermal rate calculations is the Quantum Instanton (QI) approach.^{32}^{-}^{35} Here, we only briefly sketch the key equations for QI theory and the steps for computing the relevant quantities using PIMD. The derivation of the QI rate equation begins with the formally exact quantum mechanical expression for the thermal rate constant:^{32}

(39)

where *Q _{r}*(

(40)

In the above equation, *C _{ff}* (0) is the zero time value of the flux-flux correlation function generalized to the case of two separate dividing surfaces

(41)

and Δ*H* is a particular type of energy variance

(42)

with * _{a}* and

(43)

Rewriting Eq. (40) in the form

(44)

leads to components that are easily computable using PIMD, where

(45)

(46)

(47)

The conditional average
is computed from the ensemble sampled with the *P* and *P*/2 PI slices constrained to the dividing surfaces

(48)

where the quantities within the average are defined as follows:

(49)

(50)

(51)

In the expression for *G*, *f* is the total number of degrees of freedom (i.e., *f* = 3*N* in our malonaldehyde application).

All factors required to calculate the QI rate can be computed using the PIMD facilities within Amber.^{28} For example, the joint distribution [Eq. (45)] can be computed using umbrella sampling along the reaction coordinates of the *P* and *P*/2 slices. A set of 2-dimensional (2D) biased simulations, each enhancing the sampling near a particular point of the 2D (*ξ _{P}*×

(52)

Thus, whereas TI integration with respect to mass within PIMD uses the same estimator for sampling at the dividing surface and in the RS, QI uses Eq. (52) for the sampling at the dividing surfaces and Eq. (38) for sampling in the RS, the difference being the last term in Eq. (52).

The symmetric intramolecular proton transfer reaction in malonaldehyde is depicted in Figure 1. All classical and quantum nuclear molecular dynamics simulations were performed using the DG-EVB facility within Amber 11.^{28} The DG-EVB surface was fit using the energy, gradient and Hessian information at the RS, TS, PS and off-TS geometries. The geometries at the RS, TS and PS were optimized according to the W1BD^{78}^{-}^{79} model chemistry using a development version of the Gaussian suite.^{80} The additional non-optimized off-TS geometry lies *in front* of the TS and corresponds to both OH bonds being 1.5 Å. This point was chosen to improve the surface fit to the repulsive wall in front of the TS. Being able to systematically incorporate additional ab initio points as needed for refining the surface fit is one strength of the DG-EVB approach. The gradients and Hessian tensors were computed from the resulting geometries at the CCSD/cc-pVTZ level of theory. This choice of ab initio methods is comparable to the most accurate approach described in the literature for malonaldehyde, employing frozen-core CCSD(T)/aug-cc-pV5Z single point calculations at frozen-core CCSD(T)/aug-cc-pVTZ optimized *C _{s}* and

The diabatic state Hamiltonian matrix element is approximated as a harmonic expansion about the ab initio optimized minimum (RS or PS) plus a scaled non-bonding, van der Waals, exponential-6 term from the universal force field (UFF)^{81} and Amber force field terms for angles and dihedrals

(53)

In the above equation,
and
are UFF exponential-6 parameters,
is the selected repulsive coordinate (e.g. the distance between the transferring H and the acceptor O in malonaldehyde) for *H _{NN}* and

(54)

is a constant ensuring that the ab initio energy is recovered at the DG-EVB data point, **q** = **q*** _{N}*. A set of redundant internal coordinates,

The *classical* RC *ξ*, chosen to represent the breaking and formation of a chemical bond, is defined as the difference of bond lengths *r* between the donor (D), the acceptor (A) and the transferring particle (H) positions

ξ(q) = *r*(q_{D}, q_{H})−*r*(q_{A}, q_{H})

(55)

*ξ*(**q**)<0, thus, represents the reactant state region while *ξ*(**q**) > 0 represents the product state and *ξ*(**q**) = 0 delineates the transition state (TS). Within the umbrella sampling framework,^{71} the system Hamiltonian is described by the modified potential

(56)

where **q** is the set of system coordinates, *k* is the harmonic force constant parameter and
is a biasing potential that is added to the original system potential *ε*_{0} (obtained from diagonalization of the EVB matrix) to enhance the sampling of a predetermined region of configuration space near
. The distributions of *ξ*^{(}^{n}^{)} from all of the (*n*) biased simulations are unbiased and combined using WHAM^{71}^{,}^{76} to form the PMF describing the chemical reaction of interest. All biased sampling simulations were performed in the NVT ensemble at a temperature of 300 K using a leap-frog Verlet integration time step of 0.5 fs for a minimum of 3 ns of total sampling per window.^{82} The temperature was maintained via coupling of the system to a Langevin thermostat^{83} with a collision frequency of 1 ps^{-1}.

With the exception of the modifications described below, similar simulation protocols were followed for the PIMD simulations. The RC chosen to describe the breaking and formation of a chemical bond within the PI-QTST framework is the difference of bond lengths *r* between the donor, the acceptor and transferring particle centroid coordinates

(57)

where the centroid position for particle *χ* is defined as the average of positions over the *P* path integral slices

(58)

Substituting * _{c}* for

Two procedures were utilized for estimating the primary KIE in malonaldehyde. In the direct approach, the PI-QTST^{72} rates for the intramolecular proton transfer of the hydrogen and deuterium isotopes are computed separately using Eq. (31) to form the ratio
. The centroid density component of the PI-QTST rate was computed using umbrella sampling with *k*^{(n)} = 100.0 kcal/mol-Å^{2} and
Å. The dynamical frequency factor was computed again with umbrella sampling with the RC restrained to the transition state region, i.e. *ξ*_{0} = 0.0 Å, using a harmonic force constant of *k* = 2000.0 kcal/mol-Å^{2}. A total of 14 independent biased trajectories (seven for each isotopically labeled proton transfer) are required for estimating the ratio of centroid densities, while two independent biased simulations are required for estimating the ratio of the frequency factors. In the TI with respect to mass approach,^{73}^{-}^{74} the ratio of isotopic centroid densities is estimated using Eq. (37) from a set of *λ* = {0.0, 0.2, 0.4, 0.6, 0.8, 1.0} simulations where the average _{RS} is computed from ground-state EVB sampling in the RS region and the average
is computed from umbrella sampling with the RC restrained to *ξ*_{0}=0.0 Å via a harmonic force constant, *k* = 2000.0 kcal/mol-Å^{2}. The above TI by mass, therefore, entails a total of *6 λ* × 2 trajectories/*λ* = 12 independent trajectories. Although the computational costs appear similar for both alternatives, it is expected that the statistical error for the KIE estimate based on relative rates [from Eq. (35)] will be smaller compared to the estimate based on absolute rates [from Eq. (31)]. For cases where a large number of biased trajectories are required to cover the range of RC values, due to an intrinsically high free energy barrier of the system or because of multiple dividing surfaces, the TI by mass route will be computationally less expensive.

Particularly in the case of QI^{33}^{-}^{34} calculations, TI by mass^{73}^{-}^{74} is a substantially more efficient approach for estimating KIEs. In the direct approach, the conditional averages of *f _{v}*,

Estimating the uncertainty in a rate calculation is not straightforward, especially when the quantities are derived from sampling over phase space. The frequency factor, PMF and TI by mass calculations are subject to uncertainties related to the level of convergence of phase space sampling. To provide an estimate of this type of uncertainty, we compute the various components as a parameter of MD sampling intervals. The uncertainties reported in Table 1, using this approach, give an indication of the variability of the results as a parameter of phase space averaging. Adequate MD sampling should give averages with relatively small standard deviation.

To evaluate the effectiveness of the DG-EVB method for constructing ab initio-based reactive PESs for modeling chemical dynamics, the classical TST rate, PI-QTST rate, QI rate, and KIE for the prototypical intramolecular proton transfer reaction in malonaldehyde are computed. The calculations of these observables require MD sampling of configurational space to obtain thermodynamic averages. Furthermore, since the chemistry under consideration involves a proton, nuclear quantum effects, such as zero-point motion and tunneling, are important for describing the reaction rate. These elements have been integrated into the latest release of the Amber biosimulation suite (version 11),^{28} and the results for malonaldehyde are described below.

The 2D PES for the intramolecular proton transfer reaction in malonaldehyde, employing 4 Gaussian centers (located at the RS, TS, PS and off-TS) to fit
, is depicted in Figure 2. The diabatic state *H _{NN}* is represented as a quadratic expansion about the ab initio minimum configuration, augmented by a non-bonding, van der Waals, exponential-6 term from the universal force field

Two-dimensional potential energy surface for the intramolecular proton transfer reaction in malonaldehyde employing a 4 Gaussian fit for *H*_{12} and *H*_{nn} with a UFF repulsive term and Amber angle and dihedral terms. The minima are denoted with a green dot **...**

With a well-defined PES, it is now possible to perform nuclear dynamics on this surface incorporating ab initio data. Since the ab initio energy barrier of 4.08 kcal/mol is well above thermal energy, conventional MD may not adequately sample important TS configurations. The resulting PMF (Fig. 3) obtained from ground-state dynamics indicates that malonaldehyde predominantly samples the conformational space of the reactant and product region but not so much the barrier region. Figure 3 is obtained by initiating conventional MD on the reactant minimum and collecting the statistics (indicated by the gray histogram in the background) for sampling a particular value of the RC. At an average temperature of 300K, the system has sufficient thermal energy to overcome the energy barrier. Sampling near the TS, however, remains statistically insignificant. This intrinsic *rare event* nature of the chemistry, thus, requires enhanced sampling techniques for sampling conformations within the important TS region.

Potential of mean force obtained from direct sampling on the DG-EVB ground-state surface (). The gray background depicts the histogram of RC values encountered during the MD.

Utilizing the umbrella sampling procedure^{71} described in the methodology section, a set of independent biased trajectories was used to map out the distribution of RC values over the entire range of the reactive path. Each trajectory only enhances the sampling about a predefined RC value,
. Using the WHAM^{71}^{,}^{76} procedure, the set of biased distributions (gray curves, Fig. 4) are combined to form a free energy profile ( curve, Fig. 4) spanning the entire range of the RC for the proton transfer under observation. The resulting normalized density contribution [Eq.(29)] to the TST rate is 2.04 ×10^{−3} (± 1.81 ×10^{−4}) Å^{−1}. Again, using umbrella sampling with the RC restrained to the dividing surface (*ξ*_{0} = 0.0Å), the other contribution from the dynamical frequency factor can be estimated from the gradient of the RC. With a value of 2.48 ×10^{13}(± 1.31 ×10^{9}) Å·*s*^{−1} for *dξ/dt**ξ*^{≠}, the classical TST estimate of the proton transfer rate is 2.53 ×10^{10} (± 2.25 ×10^{9})s^{-1}.

Potential of mean force obtained from umbrella sampling on the DG-EVB ground-state surface. The dashed line is a cubic spline fit through the data points (). The gray curves in the background show the histogram of RC values sampled along the **...**

While the above classical model of the proton experiences a free energy barrier of ~3.75 kcal/mol, zero-point motion and nuclear tunneling may effectively lower this reaction barrier. Figure 5 compares the quantum PMF () obtained from PIMD sampling on the DG-EVB surface to the classical PMF (dotted line) from Fig. 4. Both curves are normalized to their respective RS partition function. The free energy values at the TS (*w***(***ξ*^{≠}**)** and
), therefore, are relative to absolute energy origins as defined in the classical TST and PI-QTST rate equations. It is seen that full nuclear quantization of malonaldehyde lowers the free energy barrier by about 2 kcal/mol. This result is similar to trends observed in earlier studies of malonaldehyde by Tuckerman^{84} and by Schofield.^{85}^{-}^{86} While inclusion of nuclear quantization sufficiently lowers the free energy barrier to allow sampling of the TS, a general framework for enhancing sampling of a particular region of configurational space is desirable. Reactions in enzyme environments, for example, may encounter barriers typically on the order of ~10 kcal/mol or higher. Our implementation of umbrella sampling within the PIMD function in Amber will permit studies of rare event phenomena using quantum dynamics approaches.

From a computational efficiency perspective, it is also worthwhile to consider if nuclear quantization is necessary for the entire system or can certain degrees of freedom remain in the classical description. Figure 6 shows the PMFs as a parameter of the *level of nuclear quantization*, where the curve reproduces the classical description from Fig. 4 and the ○ curve depicts the fully quantum description. Here, all PMFs are obtained using the umbrella sampling protocol described in the Simulation Details section. The □ curve reflects quantization only of the donor, acceptor and transferring proton, while the curve shows the impact of quantization only of the transferring proton. Both of these curves are on top of each other and are slightly higher than the fully quantized system. The remaining curve (+) corresponds to quantization of only the donor and acceptor oxygens. As expected from theoretical considerations, nuclear quantization of the light particles (hydrogens) has the most significant effect on the free energy barrier height; whereas, quantization of the heavier oxygen atoms only has a minor effect compared to the classical description. For the case of malonaldehyde, quantizing the transferring proton is sufficient to capture the bulk of the quantum effects associated with our definition of the RC. This trend may not be generally transferrable to other more flexible molecules or to alternative prescriptions of the reaction path. In addition to lowering the reaction barrier, nuclear quantization also shifts the location of the reactant and product minima and changes the curvature in the transition state region. This alteration of both the barrier height and shape of the free energy surface will manifest in KIE measurements.

Potential of mean force as a parameter of the level of nuclear quantization: classical nuclei (), quantization only of the donor and acceptor oxygens (+), quantization of transferring proton (), quantization of donor and acceptor oxygens **...**

The estimates of primary KIEs in this study were obtained using both the direct approach of explicitly forming the ratio of the absolute rates and via TI with respect to mass. Table 1 shows all contributions to the chemical rates for the hydrogen and deuterium isotopes as computed from classical TST, PI-QTST and QI prescriptions. Because experimental measurements of KIEs for malonaldehyde are unavailable, the latter two approximate quantum rates serve to provide ball park estimates from well-established methods. The KIE of 4.05 (± 0.27) from PI-QTST is about one and a half times larger than the QI value of 2.41 (± 0.86). The absolute rates from PI-QTST, however, are approximately two to three times larger than the estimates from QI, due predominantly to the larger frequency factors. The contribution to the rate from the PMF is smaller in PI-QTST than in QI for both proton and deuterium transfers. Since PI-QTS and QI are approximate quantum rate methods, equivalently at the level of transition state theory, the differences in the rates are not due to dynamical recrossing effects. A clear relationship between PI-QTS and the QI method has yet to be established. Furthermore, the results indicate that classical TST is inadequate in providing KIE estimates. Here, isotopic substitution predominantly impacts the frequency factor and leaves the PMF contribution to the rate little changed, resulting in a KIE of 1.54. The effective lowering of the barrier height due to zero-point motion and nuclear tunneling are missing. Nuclear quantization via PI recovers these effects and provides the dominant contributing factor of 2.91 (± 0.19) to the PI-QTST KIE. The frequency factors obtained from PIMD are similar to those computed from classical sampling, suggesting that the centroid coordinates of the system at the dividing surface for this symmetric transfer can be represented by classical nuclear coordinates. For asymmetric reactions, this correspondence between the centroid and classical variables may not hold, and one cannot simply assume that the computationally less demanding classical MD sampling will be sufficient for estimating the frequency factor component of the PI-QTST rate.

While a direct comparison of our computed KIEs with the published results of Schofield^{85} is complicated by differences in the underlying PESs, it is worthwhile to identify the contributing factors giving rise to varying estimates of the KIE. In their previous PI-QTST study employing a molecular mechanical EVB potential for describing the proton transfer reaction,^{85} the primary KIEs range from 6.49 to 11.41 and depend on the choice of RC. The sensitivity of KIEs on the prescription of the RC is to be expected, as PI-QTST is a transition state approximation to the quantum rate. Different RCs may result in varying levels of barrier recrossings and it is this recrossing factor that is missing in the KIE calculations. Furthermore, the molecular mechanical parameters employed in that EVB formulation give rise to an energy difference of 8.75 kcal/mol between the reactant and transition state conformations. This activation energy is about two times larger than the barrier obtained from the present DG-EVB approach, leading to rates that are about 2-3 orders of magnitude smaller. Also, a higher barrier typically leads to a larger KIE, and this trend is observed between the two methods.

Figure 7(a) shows the 2D PMF associated with the joint distribution of the RC along the *P* and *P*/2 PI slices for the proton transfer reaction. The contribution to the QI rate is obtained using the dividing surface corresponding to the top of the free energy barrier in Fig. 7(b) with *ξ _{a}*=

As noted in the methods section, the computationally expensive PMF calculation can be avoided if one reformulates the ratio of isotope densities as a thermodynamic integration over mass. Figure 8 displays the average values of −*β* Φ(*λ*)/*λ* sampled during PIMD in the RS well and at the dividing surface as a parameter of *λ*. Integration over the *λ* values using Eq. (37) provides a centroid density value of 3.08 (± 0.01), which is comparable to 2.91 (± 0.19) obtained directly from the PMFs. Similarly, TI by mass for the QI calculations provide an estimate of 4.73 (± 0.02) for
compared to the value of 2.24 (± 0.06) from direct calculation. In the direct calculations, the QI PMF barrier height for proton transfer is 1.68 kcal/mol and is 2.16 kcal/mol for deuterium transfer [see Figure 7(b)]. To get a
of ~4.73, the proton transfer barrier height needs to lower to 1.23 kcal/mol (while keeping the deuterium transfer barrier at 2.16 kcal/mol) or the deuterium transfer barrier height needs to increase to 2.60 kcal/mol (while keeping the proton transfer barrier at 1.68 kcal/mol). This resolution difference of ~0.5 kcal/mol is especially challenging to achieve in the QI 2D PMFs, where errors in the reweighting of the 121 biased distributions may accumulate and impact the global unbiased distribution. The virial estimator for the average quantities in TI by mass has been shown to converge efficiently.^{74} The errors in the TI by mass estimate of the KIE, thus, are expected to be smaller than errors arising from the direct approach.

Average value of −*β**d*Φ(*λ*)/*dλ* sampled during PIMD in the RS region () and at the dividing surface
(○) as a parameter of *λ*.

Multiplying the ratio of isotopic densities by the corresponding ratio of frequency factors, the PI-QTST estimate of 4.27 (± 0.01) for the KIE is within 80% of the QI value of 5.10 (± 1.81). The error in the QI estimate, however, is much larger than that from PI-QTST because of the variability of the frequency factor estimates. While both approaches average the frequency factor over the same time span (3 ns), QI appears to require more phase space sampling to achieve a higher level of convergence. On the other hand, the phase space averages required for TI by mass converges much faster. A ns of sampling of the virial estimators for *d*Φ/*dλ* at each value of *λ* is sufficient for obtaining the ratio of isotopic densities with relative uncertainties that are less than 1%. The savings in computation as well as a much smaller estimate of the uncertainty is especially evident in QI calculations, where the two set of computationally demanding PIMD sampling can be avoided altogether by TI integration over the mass.

This work presents an application of the distributed Gaussian EVB method for constructing an ab initio-based reactive potential energy surface describing the intramolecular proton transfer reaction in malonaldehyde. Albeit a deceptively simple gas-phase system, malonaldehyde has all the salient features necessary to test the methodology development towards a description of chemical dynamics in complex systems. Extension of the approach to treat nuclear quantum effects is made possible through coupling to Path Integral methods as well as by the Path Integral Quantum TST and Quantum Instanton formalisms. Experimental measurements of rates and KIEs are important tools of chemical kinetics for elucidating the mechanism of complex chemical reactions. In situations where the KIE is important, one can either estimate this measurable by directly computing the ratio of the isotopic rates or by TI integration with respect to mass. Our results for malonaldehyde show both approaches to give similar KIEs, although TI integration with respect to mass is computationally less expensive compared to explicit construction of the PMFs. More importantly, the above functionalities of the DG-EVB approach have been made accessible to the broader community via integration within the Amber biosimulation suite.

Some needed improvements to the DG-EVB method, nevertheless, remain. The most apparent of these is our strategy for filling in the swimming hole near the transition state region by including a UFF repulsive interaction and Amber angle and dihedral terms to the diabatic state energy. Although this practical patch appears functional, a more general solution is desirable. Using a more realistic dissociative potential beyond the harmonic approximation, such as a Morse-type interaction, and using a full force field prescription for the *H _{ii}* terms will be the subjects of future development of the DG approach. It is important to emphasize that

Additionally, the current DG implementation requires the user to optimize the Gaussian *α _{K}* parameters such that the DG energy, gradient and hessian reproduce corresponding ab initio information along IRC geometries. For the case of a single DG data point, this parameter is unique. When multiple Gaussians are used to enhance the fit, the choice of

The Office of Naval Research (N00014-05-1-0457) supported this research. JLS is grateful for the computer time on the Wayne State University Grid. KFW is grateful for an allocation of computer time from the Center for High Performance Computing at the University of Utah. A portion of the computational resources for this project have been provided by the National Institutes of Health (Grant # NCRR 1 S10 RR17214-01) on the Arches Metacluster, administered by the University of Utah Center for High Performance Computing

Supporting Information **Available:** Details of the electronic structures for malonaldehyde used in the DG-EVB surface fit are provided as supporting information. This material is available free of charge via the Internet at http://pubs.acs.org.

1. Freddolino PL, Liu F, Gruebele M, Schulten K. Biophys J. 2008;94:L75. [PubMed]

2. Bader D, editor. Petascale Computing: Algorithms and Applications. Chapman & Hall/CRC Press; Boca Raton: 2008.

3. Maragakis P, Lindorff-Larsen K, Eastwood MP, Dror RO, Klepeis JL, Arkin IT, Jensen MØ, Xu H, Trbovic N, Friesner RA, Palmer AG, III, Shaw DE. J Phys Chem B. 2008;112:6155. [PMC free article] [PubMed]

4. Patel S, Brooks CL., III J Comput Chem. 2004;25:1. [PubMed]

5. Patel S, Mackerell AD, Jr, Brooks CL., III J Comput Chem. 2004;25:1504. [PubMed]

6. Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. J Comput Chem. 2004;25:1157. [PubMed]

7. Hornak V, Abel R, Okur A, Strockbine B, Roitberg A, Simmerling C. Proteins: Structure, Function, and Bioinformatics. 2006;65:712. [PubMed]

8. Xu Z, Luo HH, Tieleman DP. J Comput Chem. 2007;28:689. [PubMed]

9. Soares TA, Hünenberger PH, Kastenholz MA, Kräutler V, Lenz T, Lins RD, Oostenbrink C, van Gunsteren WF. J Comput Chem. 2005;26:725. [PubMed]

10. Hill JR, Freeman CM, Subramanian L. In: Rev Comput Chem. Lipkowitz KB, Boyd DB, editors. Wiley-VCH; New York: 2007. p. 141.

11. Balint-Kurti GG. Adv Chem Phys. 1975;30:137.

12. Eyring H. Trans Farad Soc. 1938;34:3.

13. Evans MG, Polanyi M. Trans Farad Soc. 1938;34:11.

14. Ogg RA, Polanyi M. Trans Farad Soc. 1934;31:604.

15. Evans MG, Warhurst E. Trans Farad Soc. 1938;34:614.

16. Venkatnathan A, Voth GA. J Chem Theory Comput. 2005;1:36.

17. Warshel A, Weiss RM. J Am Chem Soc. 1980;102:6218.

18. Warshel A. Computer Modeling of Chemical Reactions in Enzymes and Solutions. John Wiley & Sons, Inc.; New York: 1997.

19. Dal Peraro M, Ruggerone P, Raugei S, Gervasio FL, Carloni P. Current Opinion in Structural Biology. 2007;17:149. [PubMed]

20. Hutter J, Curioni A. Parallel Computing. 2005;31:1.

21. Carloni P, Rothlisberger U, Parrinello M. Acc Chem Res. 2002;35:455. [PubMed]

22. Bolton K, Hase WL. In: Modern Methods for Multidimensional Dynamics. Thompson DL, editor. World Scientific; Singapore: 1998. p. 143.

23. Hase WL, Song KH, Gordon MS. Computing in Science and Engineering. 2003;5:36.

24. Schlegel HB. Bull Kor Chem Soc. 2003;24:837.

25. Schlegel HB, Sonnenberg JL. Journal of Chemical Theory and Computation. 2006;2:905.

26. Sonnenberg JL, Schlegel HB. Mol Phys. 2007;105:2719.

27. Sonnenberg JL, Wong KF, Voth GA, Schlegel HB. Journal of Chemical Theory and Computation. 2009;5:949.

28. Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B. P.; Wang, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vaníček, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A.; University of California: San Francisco, 2010.

29. Feynman RP, Hibbs AR. Quantum Mechanics and Path Integrals. McGraw-Hill; New York: 1965.

30. Feynman RP. Statistical Mechanics. Benjamin; Reading, MA: 1972.

31. Kleinert H. Path Integrals in Quantum Mechanics, Statistics, and Polymer Physics. World Scientific; Singapore: 1995.

32. Miller WH. J Chem Phys. 1975;62:1899.

33. Miller WH, Zhao Y, Ceotto M, Yang S. J Chem Phys. 2003;119:1329.

34. Yamamoto T, Miller WH. J Chem Phys. 2004;120:3086. [PubMed]

35. Miller WH, Schwartz SD, Tromp JW. J Chem Phys. 1983;79:4889.

36. Wang Y, Braams BJ, Bowman JM, Carter S, Tew DP. J Chem Phys. 2008;128:224314. [PubMed]

37. Viel A, Coutinho-Neto MD, Manthe U. J Chem Phys. 2007;126:024308. [PubMed]

38. Mil'nikov GV, Yagi K, Taketsugu T, Nakamura H, Hirao K. J Chem Phys. 2004;120:5036. [PubMed]

39. Mil'nikov GV, Yagi K, Taketsugu T, Nakamura H, Hirao K. J Chem Phys. 2003;119:10.

40. Yagi K, Taketsugu T, Hirao K. J Chem Phys. 2001;115:10647.

41. Hazra A, Skone JH, Hammes-Schiffer S. J Chem Phys. 2009;130:054108. [PubMed]

42. Ben-Nun M, Martinez T. J Phys Chem A. 1999;103:6055.

43. Sewell TD, Guo Y, Thompson DL. J Chem Phys. 1995;103:8557.

44. Shida N, Barbara PF, Almolöf JE. J Chem Phys. 1989;91:4061.

45. Makri N, Miller WH. J Chem Phys. 1989;91:4026.

46. Ruf BA, Miller WH. J Chem Soc, Faraday Trans. 1988;84:1523.

47. Carrington T, Miller WH. J Chem Phys. 1986;84:4364.

48. Higashi M, Truhlar DG. Journal of Chemical Theory and Computation. 2008;4:790.

49. Albu TV, Corchado JC, Truhlar DG. J Phys Chem A. 2001;105:8465.

50. Kim Y, Cochado JC, Villa J, Xing J, Truhlar DG. J Chem Phys. 2000;112:2718.

51. Lin H, Pu J, Albu TV, Truhlar DG. J Phys Chem A. 2004;108:4112.

52. Tishchenko O, Truhlar DG. J Chem Theory Comput. 2009;5:1454.

53. Chandler D. Introduction to Modern Statistical Mechanics. Elsevier; New York: 1991.

54. Chang YT, Miller WH. J Phys Chem. 1990;94:5884.

55. Chang YT, Minichino C, Miller WH. J Chem Phys. 1992;96:4341.

56. Minichino C, Voth GA. J Phys Chem B. 1997;101:4544.

57. Saad Y, Schultz MH. SIAM J Sci Stat Comput. 1986;7:856.

58. Pulay P. Chem Phys Lett. 1980;73:393.

59. Pulay P. J Comput Chem. 1982;3:556.

60. Pulay P. In: Molecular Quantum Mechanics: Analytic Gradients and Beyond. Csaszar AG, Fogarasi G, Schaefer HF III, Szalay PG, editors. ELTE Institute of Chemistry; Budapest: 2007. p. 71.

61. Arndt M, Nairz O, Vos-Andreae J, Keller C, van der Zouw G, Zeilinger A. Nature. 1999;401:680. [PubMed]

62. Zuev PS, Sheridan RS, Albu TV, Truhlar DG, Hrovat DA, Borden WT. Science. 2003;299:867. [PubMed]

63. McMahon RJ. Science. 2003;299:833. [PubMed]

64. Chandler D, Wolynes PG. J Chem Phys. 1981;74:4078.

65. Ceperley DM. Review of Modern Physics. 1995;67:279.

66. Berne BJ, Thirumalai D. Annu Rev Phys Chem. 1986;37:401.

67. Martyna GJ, Klein ML, Tuckerman M. J Chem Phys. 1992;97:2635.

68. Hänggi P, Talkner P, Borkovec M. Review of Modern Physics. 1990;62:251.

69. Kapral R, Consta S, McWhirter L. In: Classical and Quantum Dynamics in Condensed Phase Simulations. Berne BJ, Ciccotti G, Coker DF, editors. World Scientific; Singapore: 1998.

70. Hinsen K, Roux B. J Chem Phys. 1997;106:3567.

71. Roux B. Comput Phys Commun. 1995;91:275.

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

73. Yamamoto T, Miller WH. J Chem Phys. 2005;122:044106. [PubMed]

74. Vaníček J, Miller WH. J Chem Phys. 2007;127:114309. [PubMed]

75. Cao J, Berne BJ. J Chem Phys. 1989;91:6359.

76. Kumar S, Bouzida D, Swendsen RH, Kollman PA, Rosenberg JM. J Comput Chem. 1992;13:1011.

77. Kumar S, Rosenberg JM, Bouzida D, Swendsen RH, Kollman PA. J Comput Chem. 1995;16:1339.

78. Parthiban S, Martin JML. J Chem Phys. 2001;114:6014.

79. Barnes EC, Petersson GA, Montgomery JA, Jr, Frisch MJ, Martin JML. Journal of Chemical Theory and Computation. 2009;5:2687.

80. Frisch MJ, Trucks GW, Schlegel HB, Scuseria GE, Robb MA, Cheeseman JR, Montgomery JA, Jr, Vreven T, Scalmani G, Mennucci B, Barone V, Petersson GA, Caricato M, Nakatsuji H, Hada M, Ehara M, Toyota K, Fukuda R, Hasegawa J, Ishida M, Nakajima T, Honda Y, Kitao O, Nakai H, Li X, Hratchian HP, Peralta JE, Izmaylov AF, Kudin KN, Heyd JJ, Brothers E, Staroverov V, Zheng G, Kobayashi R, Normand J, Sonnenberg JL, Ogliaro F, Bearpark M, Parandekar PV, Ferguson GA, Mayhall NJ, Iyengar SS, Tomasi J, Cossi M, Rega N, Burant JC, Millam JM, Klene M, Knox JE, 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, Chen W, Wong MW, Pople JA. Revision G.01 ed. Gaussian, Inc.; Wallingford, CT: 2007.

81. Rappé AK, Casewit CJ, Colwell KS, Goddard WA, III, Skiff WM. J Am Chem Soc. 1992;114:10024.

82. Allen MP, Tildesley DJ. Computer Simulation of Liquids. Clarendon Press; Oxford: 1987.

83. Pastor RW, Brooks BR, Szabo A. Mol Phys. 1988;65:1409.

84. Tuckerman M, Marx D. Phys Rev Lett. 2001;86:4946. [PubMed]

85. Iftimie R, Schofield J. J Chem Phys. 2001;115:5891.

86. Iftimie R, Schofield J. J Chem Phys. 2001;114:6763.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's Canada Institute for Scientific and Technical Information 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. |