Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
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/ct900579k
PMCID: PMC2992356

Proton Transfer Studied Using a Combined Ab Initio Reactive Potential Energy Surface with Quantum Path Integral Methodology


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.

1. Introduction

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, SN2 and Diels-Alder reactions. The common element among these methodologies is the description of the system within a matrix representation for the Hamiltonian


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 Hii and Hij terms. Notably, Warshel utilized a combination of molecular dynamics (MD) force fields and empirical fitting, an empirical valence bond (EVB) approach, to approximate the matrix elements for describing reactions in solutions and within enzyme environments.17-18

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 Amber28 biosimulations suite, of the DG-EVB surface generation methodology with molecular dynamics, Path Integrals,29-31 and Quantum Instanton32-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.

2. Theoretical Background

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.

2.1 Designing ab initio-based chemically reactive PESs

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


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, qTS-, the parameters A (a scalar), B (a vector) and C (a matrix) take the analytical form


When the system configuration deviates far from the transition state structure (i.e., for large Δq) and the matrix C contains one or more negative frequencies, H122 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


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


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, H122(q) is approximated as an expansion in a set of distributed Gaussians centered on a set of molecular configurations qK


where Ncfg is the number of ab initio data points used for the fitting, Ndim is the number of system coordinates, g(q, qK, i, j, αK) are the s-, p- and d-type Gaussians and BijK are the expansion coefficients. The term involving the identity matrix in Eq. (7) was accumulated into the s-type Gaussian [see Eq. (10)] to precondition the system of linear equations for efficient convergence when utilizing iterative methods. The non-standard form of the d-type Gaussian is for similar reasons. If the number of Gaussian centers, K, is equal to the number of data points where H122(q) is evaluated, Eq. (9) describes a system of linear equations


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


The F column vector stores the terms H122(qL), H122(q)/q|q=qL and 2H122(q)/q2|q=qL evaluated at the Ncfg ab initio configurations. The D matrix contains the values of the Gaussian bases, g(qL, qK, i, j, αK), [partial differential]g(q, qK, i, j, αK) /[partial differential]q[mid ]q=qL and [partial differential]2g(q, qK, i, j, αK)/[partial differential]q2[mid ]q=qL, evaluated at these same configurations qL. The B column vector contains the set of unknown expansion coefficients being solved for within this linear system of equations. Once this solution vector is obtained, we have a general (q-dependent) analytical approximation to the coupling [Eq. (9)] and corresponding derivatives [Eqs. (15) & (16)], computed as a linear combination of the distributed Gaussian basis set.

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


where ε0 is the lowest eigenvalue of the matrix and H11 and H22 are the reactant and product valence bond states. H11 and H22 can be described by a force field potential (as in the spirit of Chang-Miller54) or as a Taylor series expansion about the respective ab initio minimum.25 The key idea is to determine H122(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 ([partial differential]εΨ/[partial differential]q) and Hessian ([partial differential]2εΨ/[partial differential]q2) data in conjunction with the corresponding Hii values and derivatives. This procedure provides coupling values and derivatives for the Ncfg discrete ab initio geometries. The DG-EVB method provides a prescription for evaluating the coupling at all coordinates as a linear combination of Gaussians expanded about these ab initio configurations, which may be chosen to be distributed along the IRC (intrinsic reaction coordinate), for example. While an analytical expression of the Hij term for a general multistate N×N system Hamiltonian does not exist, one can make, nevertheless, the pair-wise approximation and estimate the couplings using the two-state expression above for unique ij pairs. In this paper, we focus only on the two-state problem.

2.2 Incorporating nuclear quantum effects

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




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


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


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

2.3 Computing observables from molecular dynamics trajectories

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


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


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/kBT, h is the Heaviside step function, ξ is the location of the dividing surface partitioning the reactant and product regions and [Xi w/ tilde] is the value of the RC that was sampled during MD described by the system Hamiltonian H(q). The dynamical frequency factor can be estimated by the velocity of a free particle along the RC direction


where [left angle bracket](...)[right angle bracket]ξ 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


in the perspective of the potential of mean force (PMF)71


where ξ0 and w(ξ0) are constants that are typically chosen to reflect initial conditions of the chemical system. For example, note that [left angle bracket]ρ(ξ)[right angle bracket] and [left angle bracket][rho with tilde](ξ)[right angle bracket] 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 [[left angle bracket][rho with tilde](ξ)[right angle bracket]]), 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


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


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 formalism72


where c/dt is the dynamical frequency factor and ρ is the density of sampled centroid RC values, defined as


where β = 1/kBT, h is the Heaviside step function, ξc is the location of the dividing surface partitioning the reactant and product regions and [Xi w/ tilde]c is the value of the RC (as a function of the centroid coordinates q(c)=1/Ps=1Pq(s)) that was sampled during MD described by the effective potential Φ(q(1), (...), q(P)) defined in Eq. (21). Similar to classical TST, the dynamical frequency factor can be estimated by the velocity of a free particle along the centroid RC direction


Furthermore, the centroid density can be recast in terms of a quantum PMF w(ξc) as


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., KIE=kPI-QTST(H)/kPI-QTST(D), 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 mass73-74


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


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


where [left angle bracket](...)[right angle bracket]RS denotes the canonical average over the reactant state region, ξc denotes the conditional average with the system constrained to the dividing surface ξc and dΦ/ can be computed from the PI virial-like estimator65,75


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Φ/ 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 ξc, 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


where Qr(T) is the reactant partition function per unit volume at temperature T, β is the inverse temperature 1/kBT and F^γ=i/[H^,h(ξγ(q)] is the flux operator where the Heaviside step function h defines the location of the dividing surface ξγ(q) = 0. Both the microcanonical density operator δ (EĤ) and the integral over E can be evaluated within the steepest descent approximation to give the following approximate expression for the QI rate constant:34


In the above equation, Cff (0) is the zero time value of the flux-flux correlation function generalized to the case of two separate dividing surfaces35


and ΔH is a particular type of energy variance


with [Delta with circumflex]a and [Delta with circumflex]b being a modified version of the Dirac delta function,


Rewriting Eq. (40) in the form


leads to components that are easily computable using PIMD, where


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


where the quantities within the average are defined as follows:


In the expression for G, f is the total number of degrees of freedom (i.e., f = 3N 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×ξP/2) configurational space is required to map out entire the QI joint distribution. Using the weighted histogram analysis method (WHAM),71,76-77 the generated biased distributions can be unbiased to form Cdd(0)/Qr on the DG-EVB ground-state surface. The remaining factors, involving the conditional average of fv, F and G, can be computed again using umbrella sampling with the P and P/2 slices constrained to the dividing surface ξ. Now, one can estimate KIEs from QI calculations by directly forming the ratio of the rates or by TI integration with respect to mass. For the latter case, the canonical average is defined over the 2D configurational space of the RS region delineated by the two dividing surfaces and the conditional average is computed with the P and P/2 slices constrained to the dividing surface ξ. Also, the modified Dirac delta functions in the conditional average of Eq. (48) introduces an additional term in the estimator for the logarithmic derivative of the delta-delta correlation function


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

3. Simulation Details

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 W1BD78-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 Cs and C2v structures.36 The present calculations yield a barrier height of 4.08 kcal/mol, nearly identical to the barrier of 4.09 kcal/mol obtained by Bowman.

Figure 1
Intramolecular proton transfer reaction in malonaldehdye.

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


In the above equation, AUFFi and BUFFi are UFF exponential-6 parameters, Δqri is the selected repulsive coordinate (e.g. the distance between the transferring H and the acceptor O in malonaldehyde) for HNN and


is a constant ensuring that the ab initio energy is recovered at the DG-EVB data point, q = qN. A set of redundant internal coordinates, q = {qr, qθ, qϕ}, comprised of bond lengths (r), angles (θ) and dihedrals (ϕ) is used in the DG-EVB method to maintain invariance of the resulting energy hypersurface under global rotations. Typically, H122(q) does not depend on the full set of redundant internal coordinates and a subspace comprised of coordinates that change within a prescribed tolerance between the reactant and product configurations has been shown to provide sufficient accuracy.26 For the malonaldehyde system, the above parameters for the angle (Kθ, θ0) and dihedral (Vn, n, δ) interactions are taken from GAFF (Generalized Amber Force Field).6 The UFF interactions were scaled by η = 0.60, and an optimized average value of 0.85 was used for all Gaussian αK parameters in this study.26

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)<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


where q is the set of system coordinates, k is the harmonic force constant parameter and Vumb(n) 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 ξ0(n). The distributions of ξ(n) from all of the (n) biased simulations are unbiased and combined using WHAM71,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 thermostat83 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


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


Substituting [Xi w/ tilde]c for ξ in Eq. (56) allows for enhanced sampling of a prescribed reaction path using the same umbrella sampling procedure as employed for the classical molecular dynamics. The quantum PMF then is generated using WHAM71,76 from the centroid densities sampled in all the biased trajectories. In contrast to the classical simulations described above, the PIMD was propagated in the normal mode representation66 using a 0.5 fs leap-frog Verlet integration time step and with the temperature maintained at 300 K via coupling to a Nosé-Hoover chain bath of size four. Each region of the reaction path in the vicinity of ξ0(n) is sampled for a minimum of 3 ns.

Two procedures were utilized for estimating the primary KIE in malonaldehyde. In the direct approach, the PI-QTST72 rates for the intramolecular proton transfer of the hydrogen and deuterium isotopes are computed separately using Eq. (31) to form the ratio kPI-QTST(H)/kPI-QTST(D). The centroid density component of the PI-QTST rate was computed using umbrella sampling with k(n) = 100.0 kcal/mol-Å2 and ξ0(n)={0.60,0.40,0.20,0.0,0.20,0.40,0.60} Å. 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 [left angle bracket](...)[right angle bracket]RS is computed from ground-state EVB sampling in the RS region and the average ξc 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 QI33-34 calculations, TI by mass73-74 is a substantially more efficient approach for estimating KIEs. In the direct approach, the conditional averages of fv, F and G [Eqs. (49)-(51)] were computed using umbrella sampling with the P and P/2 PI slices restrained to the dividing surface with a harmonic force constant of k = 2000.0 kcal/mol-Å2. The 2D joint distribution was computed using umbrella sampling with k(n) =100.0 kcal/mol-Å2 and ξ0(n)={1.00,0.80,0.60,0.40,0.20,0.0,0.20,0.40,0.60,0.80,1.00} Å for all combinations of pairs of restraints applied to the P and P/2 PI slices. Thus we have a grid of 11×11 = 121 independent biased sampling trajectories for each isotope, for a total of 242 trajectories in a single KIE estimate. Compared to the 12 trajectories required for TI by mass, the computation disparity is significant. A protocol similar to the above PI-QTST details was used for the QI TI with respect to mass procedure, with the only exception being that the conditional average in Eq. (37) is computed from umbrella sampling restraining the RC of both the P and P/2 slices to ξ0 = 0.0 Å via a harmonic force constant of k = 2000.0 kcal/mol-Å2, i.e., ξcξ(P),ξ(P/2). The relevant estimator used for the constrained sampling on the dividing surfaces is given by Eq. (52) while sampling within the RS well uses the estimator defined in Eq. (38).

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.

Table 1
Computed TST, PI-QTST and QI chemical rates and KIE for the intramolecular proton transfer reaction in malonaldehyde. The uncertainties given within parenthesis are estimated from the distribution of results as a parameter of the MD sampling intervals. ...

4. Results & Discussion

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 H122, is depicted in Figure 2. The diabatic state HNN 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 field81 and Amber angle and dihedral terms [Eq. (53)]. The UFF term was included to prevent an anomalous “swimming hole” located behind the TS region geometry.27 Alternatively, one can fill in this swimming hole by placing additional distributed Gaussians in these anomalous regions. The minima are indicated with a green dot and the TS with an orange dot. The additional configuration (indicated with a blue dot) along the plane orthogonal to the reaction path improves the curvature of the surface in this region. The goal here is to develop a robust, automated method capable of generating a multidimensional surface from ab initio information gleaned from sparse and strategically chosen geometries. For the malonaldehyde system, utilizing a quadratic expansion about the ab initio minimum and empirical FF terms is adequate for generating an accurate smooth PES for exploring chemical dynamics. The B vector of Eq. (14) was solved using the GMRES57-60 algorithm with a convergence tolerance of 1 × 10-9. The maximum error between the resulting DG-EVB and ab initio energies for geometries at the fit points was computed to be 7.49 × 10-5 kcal/mol. The corresponding maximum root mean square deviations between the DG-EVB and ab initio gradients is 7.14× 10-4 kcal/mol-Å and the maximum difference between elements of the DG-EVB and ab initio internal coordinate Hessian matrices is 5.83× 10-8 hartree/(internal coordinate)2. Our DG-EVB surface exhibits an energy barrier height of 4.08 kcal/mol, which is in excellent agreement with the best estimate of 4.09 kcal/mol published in the literature.36

Figure 2
Two-dimensional potential energy surface for the intramolecular proton transfer reaction in malonaldehyde employing a 4 Gaussian fit for H12 and Hnn 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.

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

Utilizing the umbrella sampling procedure71 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, ξ0(n). Using the WHAM71,76 procedure, the set of biased distributions (gray curves, Fig. 4) are combined to form a free energy profile ([triangle] 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 ×1013(± 1.31 ×109) Å·s−1 for [left angle bracket]dξ/dt[right angle bracket]ξ, the classical TST estimate of the proton transfer rate is 2.53 ×1010 (± 2.25 ×109)s-1.

Figure 4
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 ([triangle]). 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 ([triangle]) 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 w(ξc)), 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 Tuckerman84 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.

Figure 5
Potential of mean force obtained from direct PIMD sampling on the DG-EVB ground-state surface ([triangle]). The gray background depicts the histogram of RC values encountered during the PIMD. For comparison, the classical PMF from Fig. 4 is shown by the ...

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 [triangle] 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 [diamond] 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.

Figure 6
Potential of mean force as a parameter of the level of nuclear quantization: classical nuclei ([triangle]), quantization only of the donor and acceptor oxygens (+), quantization of transferring proton ([diamond]), 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 Schofield85 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=ξb=ξ=0.0 Å. The deuterium transfer is depicted by [triangle] and the proton transfer is depicted by ○. Note that this PMF is normalized with respect to the reactant partition function within the full 2D space and as such, one cannot simply perform a 1D biased sampling with both dividing surfaces set equal to each other. The 2D distribution requires a total of 11×11 = 121 biased sampling trajectories in order to map out the entire ξa×ξb configurational space. For estimating the KIE, the computation becomes quite expensive, totaling 242 trajectories for the malonaldehyde system.

Figure 7
Potential of mean force associated with the joint probability density Cdd(0;ξa, ξb)/Qr of Eq. (45). (a) Contours for proton transfer (0.4 kcal/mol each) as a parameter of RC values along the P and P/2 PI slices; (b) Slice through the contours ...

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 −β [partial differential]Φ(λ)/[partial differential]λ 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 Cdd(H)(0)/Cdd(D)(0) 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 Cdd(H)(0)/Cdd(D)(0) 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.

Figure 8
Average value of −β[left angle bracket]dΦ(λ)/[right angle bracket] sampled during PIMD in the RS region ([triangle]) and at the dividing surface ξc (○) 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Φ/ 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.

5. Concluding Remarks

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 Hii terms will be the subjects of future development of the DG approach. It is important to emphasize that H11 and H22 need to be higher than εΨ for the EVB approach to function satisfactorily.27 Some solutions to overcome anomalous behavior of EVB include 1) adding more ab initio data points for the fitting procedure, 2) modifying the prescription of Hii and 3) modifying the prescription of H12. The DG-EVB methods development has focused on options 1) and 2), either individually or in combination, while MCMM has explored option 3).48 In a recent paper, Truhlar and co-workers allow H12 to become imaginary (i.e. equivalent to our letting H122 become negative).52

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 αK parameters will affect the quality of the surface. Although this parameterization overhead cannot be circumvented, additional application of the DG method may provide some rule of thumb intuition as to the better choices of αK values. For example, is it necessary to assign unique values for each distributed Gaussian or will an average parameter suffice? Is it better to use αK values that provide more diffuse Gaussians or are more localized Gaussians the optimal choice? These are questions related to parameter sensitivity analysis that will be the subject of future studies. Furthermore, an extension of DG-EVB to treat chemical reactions in condensed environments requires the development of hybrid methods that fit only part of the active site region to ab initio data, as electronic structure calculations of the whole condensed phase configuration are prohibitively infeasible.

Supplementary Material



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


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. [PMC free article] [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. [PubMed]
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. [PubMed]
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. [PubMed]
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. [PubMed]
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. [PubMed]
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. [PubMed]
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.