Path-Integral Calculations of Nuclear Quantum Effects in Model

ARTICLE pubs.acs.org/JCTC

Path-Integral Calculations of Nuclear Quantum Effects in Model Systems, Small Molecules, and Enzymes via Gradient-Based Forward Corrector Algorithms Asaf Azuri, Hamutal Engel, Dvir Doron, and Dan Thomas Major* Department of Chemistry and the Lise Meitner-Minerva Center of Computational Quantum Chemistry, Bar-Ilan University, Ramat-Gan 52900, Israel

bS Supporting Information ABSTRACT: A practical approach to treat nuclear quantum mechanical (QM) effects in simulations of condensed phases, such as enzymes, is via Feynman path integral (PI) formulations. Typically, the standard primitive approximation (PA) is employed in enzymatic PI simulations. Nonetheless, these PI simulations are computationally demanding due to the large number of discretizations, or beads, required to obtain converged results. The efficiency of PI simulations may be greatly improved if higher order factorizations of the density matrix operator are employed. Herein, we compare the results of model calculations obtained employing the standard PA, the improved operator of Takahashi and Imada (TI), and several gradient-based forward corrector algorithms due to Chin (CH). The quantum partition function is computed for the harmonic oscillator, Morse, symmetric, and asymmetric double well potentials. These potentials are simple models for nuclear quantum effects, such as zero-point energy and tunneling. It is shown that a unique set of CH parameters may be employed for a variety of systems. Additionally, the nuclear QM effects of a water molecule, treated with density functional theory, are computed. Finally, we derive a practical perturbation expression for efficient computation of isotope effects in chemical systems using the staging algorithm. This new isotope effect approach is tested in conjunction with the PA, TI, and CH methods to compute the equilibrium isotope effect in the Schiff baseoxyanion ketoenol tautomerism in the cofactor pyridoxal-50 -phosphate in the enzyme alanine racemase. The study of the different factorization methods reveals that the higher-order actions converge substantially faster than the PA approach, at a moderate computational cost.

1. INTRODUCTION Enzymes are remarkably efficient catalysts evolved to perform well-defined and highly specific chemical transformations.1 Studying the nature of enzymatic rate enhancements is highly important from several aspects, including the rational design of synthetic catalysts and transition-state (TS) inhibitors. Isotope effects (IE) and particular equilibrium isotope effect (EIE) and kinetic isotope effect (KIE) are important tools in elucidating reaction mechanisms in enzymes.2 The KIE is a fundamental phenomenon measuring the sensitivity of chemical reaction rates on isotopic substitutions and provides the most direct probe to the structure of the TS of the reaction. Moreover, KIE might provide insights into tunneling in enzymes.2 EIE is an invaluable tool for insight into chemical reaction equilibrium, enzymatic binding, and hydrogen bonding.3,4 The EIE is defined as

6 is the free energy barrier. In where k is the rate constant, and ΔG¼ quantum transition-state theory (QTST),5 the exact rate constant is expressed by the QTST rate constant, kQTST, multiplied by a transmission coefficient γq:

k ¼ γq 3 kQTST

ð3Þ

where the QTST rate constant is given by kQTST ¼

1 βGðz¼6 Þ e hβQR

ð4Þ

ð1Þ

where h is Planck’s constant. In the following, we assume that γq = 1. In eq 4, G(z) is the free energy as a function of the centroid reaction 6 is the value of z[x] at the free energy maximum. coordinate z[x], z¼ Specifically, " # 6 1 Q¼ ¼ 6 Gðz Þ ¼ ln ð5Þ β ðm=2πp2 βÞ1=2

where Q is the partition function for the reactant state (RS) and product state (PS) for the light (L) and heavy (H) isotopes, and ΔGr is the reaction free energy. β = 1/kBT with kB being Boltzmann’s constant and T the temperature. Similarly, the KIE is defined as

6 is the reduced quantum phase space density at the dividing where Q¼ 6 , p = h/2π, and m is the mass. hyper-surface at z¼ The computational prediction of IE in enzymatic reactions presents a considerable challenge. First, classical statistic mechanics simulations cannot reproduce the observed KIE since

EIE ¼

KL Q PS =Q RS r r ¼ LPS LRS ¼ eβðΔGL ΔGH Þ KH QH =QH

KIE ¼

¼ 6 ¼ 6 kL eβðΔGL ΔGH Þ H k

r 2011 American Chemical Society

ð2Þ

Received: December 13, 2010 Published: April 01, 2011 1273

dx.doi.org/10.1021/ct100716c | J. Chem. Theory Comput. 2011, 7, 1273–1286

Journal of Chemical Theory and Computation they ignore nuclear quantum mechanical (QM) effects (NQE), such as zero-point energy and tunneling. Thus, computationally expensive quantum dynamics simulations are required. Second, condensed phase simulations require extensive sampling of both solute and solvent degrees of freedom to obtain converged results, and simulations are prone to statistical noise. Finally, IE depends exponentially on the free energy differences between the light and heavy isotopes, making it a very difficult observable to predict. In particular, secondary and heavy atom KIE and EIE are small in magnitude and are extremely challenging to compute from condensed phase simulations. Several simulation methods have been used to determine NQE in solution phase and enzymatic reactions. A practical approach to including these effects is via path-integral (PI) formulations which may be employed to calculate various properties of quantum or mixed quantumclassical systems.6,7 Numerous examples of PI simulations of condensed phase reactions exist.824 Additional approaches have been developed for condensed phase reactions, including the ensemble-averaged variational TS theory with multidimensional tunneling (EAVTST/MT),25,26 a wave function-based method,27,28 and model reactions.29 These methods have been applied to several enzymatic reactions with high-quality accord between the calculated and experimental KIEs. Recently, we developed a novel free energy mass-perturbation PI (PIFEP) method,13,14 which has been successfully applied to numerous model and enzymatic reactions.3034 In particular, a combined PIFEP and EA-VTST/ MT study has recently identified enhanced tunneling in the enzyme nitroalkane oxidase compared to the analogues uncatalyzed reaction.34 Recently, additional approaches for computation of IE have been developed and applied to various chemical systems.35,36 The modeling of IE, however, is computationally extremely demanding, and it is important to develop enhanced methods. In PI simulations of enzymes the standard primitive approximation (PA) is typically employed. The PA is based on a primitive factorization of the canonical density operator and when combined with efficient sampling schemes, such as staging,37 bisection,38 or normal mode transformation,39 yields satisfactory results. However, the application of these PA PI simulations to condensed phase reactions, employing fully QM or hybrid QM/molecular mechanical (MM) potential energy functions, is computationally demanding. The efficiency of PI simulations may be greatly improved if more accurate factorization schemes are employed. A higher-order PI approach was devised by Takahashi and Imada (TI) which may greatly enhance the efficiency of PI simulations.40,41 Suzuki has formulated a higher-order composite factorization scheme which features an additional squared force term.42,43 This strategy has been employed in studies targeting time-dependent classical dynamics,44 real-time propagator,45,46 solution of the FokkerPlanck equation,47 and in PI simulations.4856 Herein, we compare the results of model calculations obtained employing the standard PA, the TI,40 and a novel higher-order factorization based on the symplectic algorithms developed by Chin.46,57 These Chin-based factorizations have recently been employed in the study of quantum liquids.58 In the current study, the quantum partition function is computed for several model one-dimensional systems, including the harmonic oscillator, Morse potential, and symmetric and asymmetric double wells. The computations employ two complementary methods free of sampling noise: the numerical matrix multiplication5961 and the

ARTICLE

Scheme 1. KetoamineEnolimine Tautomerism in a Model PLP System

Scheme 2.

D-Ala

PLP Employed in AlaR

direct matrix diagonalization methods.62 Results emerging from this study on the different higher-order factorization methods reveal that the higher-order actions converge substantially faster than the PA approach, at a moderate computational cost. Moreover, we obtain a unique parametrization for the Chin factorization (CH) which is equally applicable to all the potentials employed herein. As a test case, we employ these higher-order Feynman PI formulations of the density matrix in conjunction with density functional theory (DFT) calculations to estimate the QM correction to vibrational free energy in a water molecule. Finally, we derive a perturbation expression for efficient computation of IEs in chemical systems using the staging algorithm. This new IE approach is tested in conjunction with the PA, TI, and CH methods to compute the EIE in the cofactor pyridoxal-50 -phosphate (PLP) Schiff base ketoenol tautomerism (Schemes 1 and 2) in the enzyme alanine racemase (AlaR).

2. THEORY 2.1. Basic Formal Expressions. The PI strategy is particularly suited for computing the quantum partition function Q in condensed phase systems since it may be obtained by applying classical simulation techniques.63 The partition function is defined as the trace of the canonical density matrix:

Z Q ¼ TrðFÞ ¼

Z dxFðx, x; βÞ ¼

dxÆxjeβH jxæ

ð6Þ

where x denotes the position of a particle in one dimension and extension to N dimensions is straightforward. To express the partition function as a Feynmann PI, we write the density matrix operator as a product of P exponents, each representing a time slice of length τ = β/P: Z Q ¼ 1274

dxÆxjeτH eτH 3 3 3 eτH jxæ

ð7Þ

dx.doi.org/10.1021/ct100716c |J. Chem. Theory Comput. 2011, 7, 1273–1286

Journal of Chemical Theory and Computation and insert a complete set of P 1 eigenstates, Z Q ¼ Z Z ¼ Z

dx1 Æx1 jeτH

dx2 jx2 æÆx2 jeτH

dxi|xiæÆxi| = 1:

Z dx3 jx3 æÆx3 j 3 3 3

dxP jxP æÆxP jeτH jxP þ 1 æ dx1 dx2 3 3 3 dxP Æx1 jeτH jx2 æÆx2 jeτH jx3 æ 3 3 3 ÆxP jeτH jxP þ 1 æ dx1 dx2 3 3 3 dxP Fðx1 , x2 ; τÞFðx2 , x3 ; τÞ 3 3 3 FðxP , xP þ 1 ; τÞ

¼ Z ¼

Z

R

ARTICLE

dx1 3 3 3 dxP

P Y

Fðxi , xi þ 1 ; τÞ

ð8Þ

i¼1

where x1 = xPþ1. In the limit Pf¥ and τf0, one can use the semiclassical PA: Fðxi , xi þ 1 ; τÞ = FPA ðxi , xi þ 1 ; τÞ ¼ FT ðxi , xi þ 1 ; τÞFV ðxi ; τÞ

ð9Þ

where FT is the kinetic energy (T), i.e., free particle term: m 2 FT ðxi , xi þ 1 ; τÞ ¼ Ω 3 exp τ 2 2 ðxi xi þ 1 Þ ð10Þ 2τ p where Ω = (m/2πτp2)1/2 and m is the mass, while FV is the potential energy (V) term: FV ðxi ; τÞ ¼ exp½τV ðxi Þ

The need for a correction term arises due to the fact that the kinetic and potential energy operators do not commute. Thus, in order to add the TI correction, all that is required is to compute the gradient of the potential. Here we suggest employing a family of higher-order factorization methods based on the symplectic algorithms developed by Chin.46,57 We start with the expression suggested by Chin (eq 29 in ref 46): eτðT þ V Þ = et3 τT ev3 τV ða3 τÞ et2 τT ev2 τWða2 τÞ et1 τT ev1 τV ða1 τÞ et0 τT

eτðT þ V Þ = ev3 τV ða3 τÞ et2 τT ev2 τWða2 τÞ et1 τT ev1 τV ða1 τÞ e2t0 τT ð17Þ Substituting for W yields: expfτðT þ V Þg = expfv3 τV ða3 τÞg expft2 τTg exp v2 τðV ða2 τÞ þ ðu0 =v2 Þτ2 ½V ða2 τÞ, ½T, V ða2 τÞÞ expft1 τTg expfv1 τV ða1 τÞg expf2t0 τTg ð18Þ If the computation of the commuter in eq 18 is not the bottleneck of the calculation (e.g., QM/MM simulations), it is advantageous to distribute the commuter more evenly over the three V. Thus, we multiply the central commuter term by a factor of 1 λ and add λ/2 times the commuter term to each potential operator on each side, as suggested by Chin,44,46,65 obtaining: expfτðT þ V Þg = exp v3 τ V ða3 τÞ þ ðλu0 =2v3 Þτ2 ½V ða3 τÞ, ½T, V ða3 τÞ expft2 τTg exp v2 τ V ða2 τÞ þ ðð1 λÞu0 =v2 Þτ2 ½V ða2 τÞ,

ð12Þ

where O(τ3) is the big O notation describing the convergence as a function of τ. Summation of the exponential in computing a property (e.g., partition function) P times, yields an error order of O(τ2) (i.e., the error in the partition function computed using PA decreases quadratically as P increases or β decreases). In the higher-order action due to Takahashi and Imada (TI),40 the following operator decomposition is employed:

ð16Þ

where W is an effective potential given by W = V þ (u0/ v2)(τ2[V,[T,V]]) and al, tl, and vl are positive coefficients which will be defined explicitly below. Setting t3 = t0 in eq 16 and redistributing the kinetic energy term at t0:

ð11Þ

where V (xi) is the potential at time slice i. The above expression (eq 9) is correct in the Pf¥ limit due to the Trotter formula eβ(TþV) = limPf¥(eτTeτV)P.64 The above quantum system is isomorphic to a classical system of ring polymers where each bead, i, in the polymer interacts with its neighbor, i ( 1, via a harmonic potential (eq 10) and experiences only a fraction, 1/P, of the full potential V (eq 11). In order to obtain improved PI methods, various operator decomposition approaches may be employed. First, we note that the PA algorithm may be derived from the general operator splitting: eτðT þ V Þ ¼ eτV =2 eτT eτV =2 þ Oðτ3 Þ

where WTI(xi) is the effective one-dimensional TI potential: p2 τ2 jFðxi Þj2 where WTI ðxi Þ ¼ V ðxi Þ þ 24m DV ðxi Þ Fðxi Þ ¼ ð15Þ Dxi

½T, V ða2 τÞÞg expft1 τTg

exp v1 τ V ða1 τÞ þ ðλu0 =2v1 Þτ2 ½V ða1 τÞ, ½T, V ða1 τÞ expf2t0 τTg

ð19Þ

Setting a1 = t0, a2 = 1/2, a3 = 1 t0, and v3 = v1 we get expfτðT þ V Þg = exp τ v1 V ðð1 t0 ÞτÞ þ ðλu0 =2Þτ2 ½V ðð1 t0 ÞτÞ,

expfτðT þ V Þg = expfτTg expfτðV þ τ2 ½V , ½T, V =24Þg

ð13Þ

½T, V ðð1 t0 ÞτÞÞg expft2 τTg exp τ v2 V ðτ=2Þ þ ð1 λÞu0 τ2 ½V ðτ=2Þ, ½T, V ðτ=2Þ

This expression converges as O(τ5) with respect to the diagonal terms (e.g., the partition function may be computed with fourthorder accuracy). This yields the density matrix elements: m FTI ðxi , xi þ 1 ; τÞ ¼ Ω 3 exp τ 2 2 ðxi xi þ 1 Þ2 τWTI ðxi Þ 2τ p ð14Þ

expft1 τTg

exp τ v1 V ðt0 τÞ þ ðλu0 =2Þτ2 ½V ðt0 τÞ, ½T, V ðt0 τÞ expf2t0 τTg 1275

ð20Þ


Journal of Chemical Theory and Computation

ARTICLE

Substituting for the commuter [V,[T,V]] = (p2/m)|F|2, where F is the gradient as defined above in eq 15, yields eτðT þ V Þ ¼ ev1 τWi et2 τT ev2 τWj et1 τT ev1 τWk e2t0 τT

integral which may be evaluated directly or simulated using Monte Carlo, MC, methods). This may be seen by inspecting eq 23. This CH PI approach with a gradient-based forward correction converges as τ6 in favorable cases, such as a harmonic potential, compared with the τ4 convergence of TI and τ2 for PA. 2.2. Condensed Phase Expressions. In condensed phase simulations it is useful to compute the QM effects as a correction to the classical mechanics (CM) results. Thus, we write the ratio between the classical and quantum partition functions:9,10 Z dxFQM ðx, x; βÞ Q QM ¼Z ð26Þ Q CM dxFCM ðx, x; βÞ

ð21Þ

This expression represents a family of algorithms with fourth-order convergence, which may be modified by changing the parameters u0, v1, v2, t0, t1, and t2.46 Interestingly, O(τ6) convergence may be achieved by an optimal choice of factorization parameters, due to cancellation of higher-order error terms. Based on this CH, the corresponding density matrix then becomes

Z 1 1=2 dxj dxk FCH ðxi , xi þ 1 ; τÞ ¼ Ω3 3 3 2t12 t0

8

9 1 1 1 > 2 2 2 > < τ m = ðx x Þ þ ðx x Þ þ ðx x Þ i j j k k i þ 1 t1 2t0 2τ2 p2 t1 exp > > : τðWðxi Þ þ Wðxj Þ þ Wðxk ÞÞ ;

Here the QM density matrix may be described by PA, TI, or CH, as described above, while the CM density matrix may be written as an analogue of the PA, TI, and CH approaches, respectively. In general, we may write the high-temperature density matrices:

ð22Þ

FQM ðxi , xi þ 1 ; τÞ ¼ FT ðxi , xi þ 1 ; τÞFΜ V ðxi ; τÞ

ð27Þ

where i, j, and k correspond to time slices (1 t0)τ, τ/2, and t0τ, respectively, and W(xi/j/k) are generalized effective one-dimensional TI-like potentials at time slices i, j, and k:

FCM ðxi , xi þ 1 ; τÞ ¼ FT ðxi , xi þ 1 ; τÞFPA V ðxc ; τÞ

ð28Þ

p 2 u0 λ jFðxi Þj2 Wðxi Þ ¼ v1 V ðxi Þ þ τ2 2m p2 u0 ð1 λÞ jFðxj Þj2 Wðxj Þ ¼ v2 V ðxj Þ þ τ2 m p2 u0 λ jFðxk Þj2 Wðxk Þ ¼ v1 V ðxk Þ þ τ2 2m

where Μ represents the PA, TI, or CH methods and xc is the classical coordinate which coincides with the centroid, x, which in discrete representation is defined as x ¼ P1∑Pi¼ 1 xi . Further we may write Z

ð23Þ

Q QM ¼Z Q CM

¼Z

Z dx1 3 3 3 dxP δðxc xÞ

Z

Z dxc

¼Z

dx1 3 3 3 dxP δðxc xÞ

Z dxc


P Y i¼1 P Y

FQM ðxi , xi þ 1 ; τÞ FCM ðxi , xi þ 1 ; τÞ

i¼1 P Y i¼1 P Y i¼1

FT ðxi , xi þ 1 ; τÞFΜ V ðxi ; τÞ FT ðxi , xi þ 1 ; τÞFPA V ðxc ; τÞ

ð29Þ where FT(xi,xiþ1;τ) and FΜ V (xi;τ) have been defined above. The delta function, δ(xc x), imposes the centroid constraint on the beads, assuring that the centroid coincides with the classical position. The classical analogue of the quantum potential energy density matrix is obtained in the limit P = 1 and is defined as FPA V (xc;τ) = exp[τV(xc)]. Employing either the PA, TI, or CH potentials, the following useful expression may be derived

ð25Þ

Z dxc FPA V ðxc ; τÞ

Q QM ¼ Q CM


dxc

This latter expression may be a useful starting point for reducing the fourth-order error. Herein λ will be limited to values between 0 and 1 to yield a forward algorithm (i.e., a negative exponent in eq 22 yields an expression with a bounded Z

Z dxc

and λ is a function of t0 yielding an algorithm correctable to sixth order for the harmonic oscillator:65 1 þ 6t0 f3 þ 4t0 ½6 þ t0 ð23 þ 24t0 Þg 5½1 12t0 ð1 2t0 Þ2 ½1 6t0 ð1 þ 2t0 4t02 Þ

dxFCM ðx, x; βÞ

Z

Here u0, v1, v2, t1, and t2 are parameters to be optimized via t0:

1 1 1 p ffiffi ffi 1 0 e t0 e ; t1 ¼ t2 ¼ t0 2 2 3 1 v1 ¼ ; v2 ¼ 1 2v1 ; 6ð1 2t0 Þ2 " # 1 1 1 1 u0 ¼ þ ð24Þ 12 1 2t0 6ð1 2t0 Þ3

λ¼

dxFQM ðx, x; βÞ

P P Y dx1 3 3 3 dxP δðxc xÞ FT ðxi , xi þ 1 ; τÞðexpðτð ðW Μ ðxi Þ V ðxc ÞÞÞÞÞ i¼1 i¼1 Z Z P Y dxc FPA dx1 3 3 3 dxP δðxc xÞ FT ðxi , xi þ 1 ; τÞ V ðxc ; τÞ

∑

i¼1

¼

Æexpðτð

P

ðW Μ ðxi Þ V ðxc ÞÞÞÞæT , x ∑ i¼1

where WM is the effective potential according to PA, TI, or CH.

c

ð30Þ V ðxc Þ

In eq 30 the internal bracket, Æ 3 3 3 æT,xc, is an average over the 1276



ARTICLE

free-particle distribution which is constrained to the centroid (classical) position, while the external average, Æ 3 3 3 æV(xc), is over the classical (centroid) potential. In this formulation, which is an extension of the original quantized classical path methods,9,10 the sampling of the classical centroid coordinate and the quantum PI coordinate may be performed separately. Enhanced sampling may be obtained by using the MC staging algorithm37 in conjunction with the expression in eq 30. In the case of M = CH, a symmetrized version of the potential must be employed in conjunction with the standard staging algorithm. 2.3. Perturbation Expression for Accurate Isotope Effects. In principle, one can carry out separate centroid path integral (PI) simulations to make QM corrections to the classical potential of mean force for different isotopes. Then, one can use the free energies for different isotopic reactions to compute the corresponding IEs. However, the statistical errors associated with these separate calculations are at least one order of magnitude greater than the free-energy difference for different isotopic reactions—an error too large to be useful for computing IEs. Thus, a sampling scheme which avoids separate sampling for different isotopes is of great importance. Here we present such a scheme for the staging algorithm. Assuming we want to sample P 1 beads using the staging algorithm, {x2, ..., xP}, between end-points x1 and xPþ1. We define x1 = xPþ1 = 0 and Λm=(2πΩ2)1/2. Stage 1: rffiffiffiffiffiffiffiffiffiffiffi xP þ 1 þ x1 ðP 1Þ P1 þ Λm η1 x2 ¼ P P xP þ 1 þ x1 ðP 1Þ þ Λm θ1 ¼ Λm θ1 ð31Þ P pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where θ1 ¼ η1 ðP 1Þ=P and η1 is a random number with normal distribution, zero mean, and unit variance. Stage 2: rffiffiffiffiffiffiffiffiffiffiffi xP þ 1 þ x2 ðP 2Þ P2 þ Λm η2 x3 ¼ P1 P1 ¼

P2 ¼ Λm θ1 þ Λm θ2 P1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where θ2 ¼ η2 ðP 2Þ=ðP 1Þ Stage 3: xP þ 1 þ x3 ðP 3Þ þ Λm η3 x4 ¼ P2

We may write xk in a more compact form: xk ¼ Λ m

We thus obtain the following identity for the free particle density matrices of the two isotopes: mL 2 ΩL 3 exp τ 2 2 ðxi, L xi þ 1, L Þ 2τ p mH 2 ¼ ΩH 3 exp τ 2 2 ðxi, H xi þ 1, H Þ ð37Þ 2τ p This is in accord with our previous work employing the bisection sampling algorithm.13 We may then write the ratio between the QM partition functions for different isotopes (i.e., IE) as IE ¼

QLQM QHQM Z Z P Y dxc dx1, L :::dxP, L δðxc xÞ FLT ðxi, L , xi þ 1, L ; τÞFΜ V ðxi, L ; τÞ

¼Z

dx1, H :::dxP, H δðxc xÞ

Z

Z

¼

Z

dx1, L :::dxP, L δðxc xÞ

Z dxc

¼Z

and and

i¼1

Z dxc

Z

ð34Þ

P Y


dxc

rffiffiffiffiffiffiffiffiffiffiffi P3 P2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where θk 1 ¼ ηk 1 ðP k þ 1Þ=ðP k þ 2Þ k1 xk 1 ¼pΛffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m ∑i ¼ 2 θi ðP k þ 1ffiÞ=ðP i þ 1Þ θi ¼ ηi ðP iÞ=ðP i þ 1Þ.

i¼1

Z dxc

RP

xk 1 ðP k þ 1Þ ¼ þ Λm θk 1 Pkþ2

ð35Þ

Thus, we see that the final bead distribution is independent of the initial position and may be written exclusively as a function of mass and random distribution numbers. In practice, we implemented the staging algorithm employing eqs 3134. Considering a reaction where the light atom of mass mL is replaced by a heavier isotope of mass mH, we use exactly the same sequence of random numbers, that is, displacement numbers {θi}, to generate the staging PI distribution for both isotopes. Thus, the resulting coordinates of these two bead distributions differ only by the ratio of the corresponding masses, assuming we use an identical random number series for the two isotopes: rffiffiffiffiffiffiffi xmk L mH ¼R ð36Þ ¼ xmk H mL

ð32Þ

P3 P3 þ Λm θ2 þ Λm θ3 ¼ Λm θ1 ð33Þ P1 P2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where θ3 ¼ η3 ðP 3Þ=ðP 2Þ In general, we may write for stage k 1: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xP þ 1 þ xk 1 ðP k þ 1Þ Pkþ1 þ Λm ηk 1 xk ¼ Pkþ2 Pkþ2

Pkþ1

k

∑ θi i¼2 P i þ 1


Z dxc


Μ FH T ðxi, H , xi þ 1, H ; τÞFV ðxi, H ; τÞ

P Y

i¼1 P Y i¼1

P Y

i¼1 P Y i¼1

FLT ðxi, L , xi þ 1, L ; τÞFΜ V ðxi, L ; τÞ

Μ FH T ðRxi, L , Rxi þ 1, L ; τÞFV ðRxi, L ; τÞ

FLT ðxi, L , xi þ 1, L ; τÞFΜ V ðxi, L ; τÞ

FLT ðxi, L , xi þ 1, L ; τÞFΜ V ðRxi, L ; τÞ

ð38Þ where we have used the substitution dxi,H = Rdxi,L. This expression may then be employed to compute IEs (e.g., EIE = IEPS/IERS).

3. COMPUTATIONAL DETAILS For the model one-dimensional systems, the numerical integration is performed with the iterative scheme for numerical matrix multiplication (NMM) as it avoids the numerical noise inherent to sampling methods, such as MC integration.60,61 Additionally, the NMM method allows rapid analysis of the convergence as a function of number of beads. The partition function is obtained by computing the trace of the density matrix 1277



ARTICLE

(eq 6). The high-temperature density matrix is computed with the PA, TI, and CH approaches. To obtain further insight into the properties of the solutions, we employ density matrix diagonalization (DMD) which gives the eigenvalues and eigenvectors of the density matrix as well as the trace.62 The partition function is computed for various well-studied potentials relevant for modeling chemical reactions: harmonic oscillator (HO), Morse oscillator (MO), symmetric double well (SDW), and asymmetric double well (ADW) which possess quantum behavior such as zero point energy and tunneling in the temperature range of 100500 K. The HO is given by VHO ðxÞ ¼ ax2

ð39Þ

where we have employed the mass of a hydrogen atom and a = 309 kcal/mol 3 Å2. The MO is given by VMO ðxÞ ¼ De ½1 eRðx x0 Þ 2

ð40Þ

where De = 136.3 kcal/mol, R = 2.2112 Å1, and x0 = 0.9166 Å. The values were chosen to resemble those of the HF molecule, and the reduced mass of HF was employed.66 The SDW is given by VSDW ðxÞ ¼ ax4 þ bx2 þ c

ð41Þ

where we have employed a mass of 1224.259 me (where me is the mass of an electron), a = 0.01, b = 0.01, and c = 0.0025 au. The ADW is given by VADW ðxÞ ¼ ax4 þ bx2 þ cx þ d

ð42Þ

where we have employed a mass of 1224.259 me, a = 0.01, and b = 0.02, c = 0.005, and d = 0.015 au. For all of the above potentials, the optimal CH value of λ (eq 25) was √ obtained by varying the parameter t0 in the range 0 to (1 1/ 3)/2. Specifically, 10 values at equal intervals were chosen: 0.0211, 0.0422, 0.0633, 0.0844, 0.1055, 0.1266, 0.1477, 0.1688, 0.1899, and 0.2110. Of these values, t0 = 0.1899 does not fall in the forward range, as it yields a positive exponent in eq 22. We also attempted to optimize λ and t0 separately.58 This was done by initially setting λ = 0 and finding the optimal t0 value. Subsequently an optimal λ value for this t0 was sought after. We found that for the potentials examined here we obtain very similar results with both approaches, and we prefer the simplicity of using eq 25. For simulations employing eq 30, we employed the CHARMM program.67 Previously, we have implemented the PA method within CHARMM.1114 In this work we also implement the TI and CH approaches together with the staging algorithm. Calculations for the water molecule employed the B3LYP functional68,69 with the 6-31þG(d,p) basis set.70 In this case, CHARMM was combined with the Gamess-UK electronic structure program.71 Simulations on the enzyme AlaR employed a hybrid QM/MM potential, where the QM part was described by a specific reaction parameter version of the semiempirical AM1 Hamiltonian.31,32 Details of the system setup and the classical molecular dynamics simulations have been published previously.31,32 In the current study we employed eqs 1 and 38 to compute IEs at a temperature of 298 K, as implemented in a development version of CHARMM. In all simulations the bead sampling was performed by simultaneously moving all beads at each PI step using the staging or mass perturbation staging

algorithms. The number of classical configurations employed was 5200, while 10 MC PI steps were performed at each classical configuration. The programs employed for all calculations on model systems are written using the Fortran programming language on a Linux platform with Intel compilers. All mathematical derivations are verified using the Maple 12 software suite (Waterloo Maple Inc.).

4. RESULTS To demonstrate the performance of the higher-order method, we apply the CH algorithm to compute the partition function (eq 6) for a number of well-studied potentials which model key features of chemical reactions. In particular, we compute the partition function for the HO, MO, SDW, and ADW. The results of CH are compared with those obtained with the PA and TI approaches. All potentials were studied at temperatures of 100, 200, 300, 400, and 500 K. A series of 10 t0 parameters were tested for the CH algorithm with all the potentials at all the temperatures studied. For the sake of brevity, only the results at 100, 300, and 500 K are presented here. The results are compared with the exact results for the HO and MO as well as for SDW and ADW.72 In order to assess the performance of the CH methods in simulations, we tested the PA, TI, and CH methods on a water molecule as well as in the enzyme AlaR. The PA, TI, and CH methods are employed for a water molecule treated with DFT and a hybrid QM/MM potential for the enzyme AlaR. Below the following notation will be employed: P refers to the number of discrete points in the PI, while k-level is defined by the integer k as P = 2k. 4.1. Harmonic Oscillator (HO). The HO serves as a simple model for a chemical bond. Key results are shown in Figure 1 and Table 1. In Figure 1 the convergence of the CH algorithm with different t0 values is presented at temperatures 100, 300, and 500 K. The optimal t0 values are 0.1055 and 0.1266 for all three temperatures, with errors increasingly greater when deviating from these optimal numbers. Indeed, these values are near the midpoint of the range given in eq 24. At 100 K the parameter t0 = 0.1899 is clearly an outlier, which is symptomatic of a positive sign in the exponent in eq 21. In Table 1 the partition function is displayed at temperatures 100, 300, and 500 K for the PA, TI, and CH. Inspection of the results at T = 100 K shows that using PA the partition function does not converge to within 1% of the exact partition function value of 3.97 109 with a k-level of up to 6. Indeed, PA converges only at a k-level of 9, corresponding to 29 = 512 integrals (results not shown in table). Using TI reduces this to k = 6, corresponding to 26 = 64 integrals, whereas with CH the desired accuracy is reached with k = 4 with the optimal t0, requiring 3 3 24 = 48 integrals. At T = 300 K the exact partition function value is 1.58 103. PA converges with k = 7, TI requires k = 4, and CH converges with k = 2. At T = 500 K the exact value of the partition function is 2.09 102. To reach convergence PA requires k = 5, TI requires k = 3, whereas CH needs k = 1. Thus, at all temperatures, the CH method yields a 25% enhancement in performance compared to TI. Moreover, CH reaches the performance of PA at less than 10% of the cost at T = 100 and 300 K, and at 500 K, it reaches the performance of PA at 20% of the cost, where we assume for simplicity that the computation of gradients does not significantly increase the computational cost. In Table 2 we present the computed IE on the partition functions of hydrogen and deuterium (QH/QD). The exact value 1278



ARTICLE

is 3.485 103 at a temperature of 100 K. It is clear from inspection of these results that PA has not yet converged (to within 1% of the exact value) at a k-level of 6, whereas TI converges with a k-level of 6. On the other hand CH reaches convergence with a k-level of 4 and 5 with use of the optimal t0 value of 0.1266 or the symmetric t0 value of 1/6, respectively. As expected, the IE converges slightly faster than the absolute partition functions. 4.2. Morse Oscillator (MO). The Morse potential is employed as a simple model for the chemical bond including anharmonicity.73 Main results are shown in Figure 2 and Table 3. In Figure 2 the convergence of the CH algorithm with different t0 values is presented at temperatures 100, 300, and 500 K. The optimal t0 values are 0.1055 and 0.1266 for all three tempera-

tures. Again, the value t0 = 0.1899 yields slightly greater errors than the other parameter values. In Table 3 the partition function is displayed at temperatures 100, 300, and 500 K for the PA, TI, and CH. Inspection of the results at T = 100 K shows that using PA the partition function does not converge to within 1% of the exact partition function value of 2.98 1013 at a k-level of 6. Rather a k-level of 10 is required using PA, corresponding to 1024 integrals (results not shown in table). Using TI reduces this to k = 7, corresponding to 128 integrals (results not shown in table), whereas with CH the desired accuracy is reached with k = 5 with the optimal t0, requiring 96 integrals. At T = 300 K the exact partition function value is 6.68 105. PA converges with k = 8, TI requires k = 5, and CH converges with k = 3. At T = 500 K the exact value of the partition function is 3.12 103. To reach convergence PA requires k = 6, TI requires k = 4, whereas CH needs k = 2. Thus, at all temperatures, the CH method yields a 25% enhancement in performance compared to TI. Moreover, CH reaches the performance of PA at less than 10% of the cost at T = 100 and 300 K, and at T = 500 K, it reaches the performance of PA at 20% of the cost. 4.3. Symmetric Double Well (SDW). The double well potential is a simple model for chemical reactions and hydrogen bonding, such as proton transfer or hydrogen networking in ice.74 Key results for this potential are shown in Figure 3 and Table 4. In Figure 3 the convergence of the CH algorithm with different t0 values is presented at temperatures 100, 300, and 500 K. The optimal t0 values are 0.1055 and 0.1266 for all three temperatures. At a k-level of 1, t0 = 0.1899 failed to converge, something we ascribe to the aforementioned sign problem. In Table 4 the partition function is displayed at temperatures 100, 300, and 500 K for the PA, TI, and CH. Inspection of the results Table 2. Isotope Effect (H/D) for the HO Calculated by the PA, TI, and CH Algorithms at Various k-Levels at T = 100 K QH/QD (T = 100 K) CHb

PA

1

5.031 101

2.584 101

2.969 102

2.977 102

2.609 10

1

2

6.384 10

3

6.484 103

8.363 10

2

3.557 10

3

3.688 103

1.938 10

2

3.414 10

3

3.492 103

3

3.486 103 3.486 103

2 3 4 5 6

Figure 1. Partition functions for the HO calculated by the CH algorithm at T = 100, 300, and 500 K, with varying values of the parameter t0.

CHa

log2 P

TI

3

6.558 10 4.176 103

7.875 10

2

1.561 10

3

4.925 10

3

3.619 10 3.495 103

3.460 10 3.478 103

a The CH method employed t0 = 1/6. b The CH method employed t0 = 0.1266.

Table 1. Partition Functions for the HO Calculated by the PA, TI, and CH Algorithms at Various k-levels at Temperatures T = 100, 300, and 500 K Q (T = 100 K)

Q (T = 500 K) a

CHa

log2 P

PA

1

2.658 103

8.303 105

1.317 107

2.297 102

5.330 103

1.701 103

5.935 102

2.811 102

2.099 102

4

6

9

3

3

3

2

2

2.090 102

2

2.089 102

2 3

a

Q (T = 300 K) a

TI

1.096 10

6

2.423 10

CH

1.470 10

8

3.947 10

PA

8.955 10

9

4.244 10

TI

6.597 10

3

2.743 10

CH

2.194 10

3

1.653 10

PA

1.589 10

3

1.584 10

TI

3.236 10

2

2.396 10

2.189 10 2.098 10

4 5

8

7.951 10 1.100 108

9

6.271 10 4.164 109

9

3.982 10 3.974 109

3

1.864 10 1.653 103

3

1.589 10 1.584 103

3

1.584 10 1.584 103

2

2.168 10 2.109 102

2

2.090 10 2.089 102

2.089 102 2.089 102

6

5.273 109

3.987 109

3.973 109

1.601 103

1.584 103

1.584 103

2.094 102

2.089 102

2.089 102

The CH algorithm employed t0 = 0.1266. 1279



ARTICLE

at T = 100 K shows that using PA the partition function converges to within 1% of the exact partition function value of 2.57 103 at a k-level of 7, corresponding to P = 128 integrals (results not shown). Using TI reduces this to k = 5, corresponding to 32 integrals, whereas with CH the desired accuracy is reached with k = 3 with the optimal t0, requiring 24 integrals. At T = 300 K the exact partition function value is 1.75 101. PA converges with k = 5, TI requires k = 3, and CH converges with k = 1. At T = 500K the exact value of the partition function is 4.59 101. To reach convergence PA requires k = 4, TI requires k = 2, whereas CH needs k = 0, the latter corresponding to three integrals. Thus, at T = 100 and 300 K, the CH method yields a 25% enhancement in performance compared to TI, while at T = 500 K TI is somewhat more efficient. Moreover, CH reaches the

performance of PA at ca. 20% of the cost at T = 100 and 300 K, while at T = 500 K PA is nearly three times as costly as CH. 4.4. Asymmetric Double Well (ADW). Principal results are shown in Figure 4 and Table 5. In Figure 4 the convergence of the CH algorithm with different t0 values is presented at temperatures 100, 300, and 500 K. The optimal t0 values are 0.1055 and 0.1266 for all three temperatures. In Table 5 the partition function is displayed at temperatures 100, 300, and 500 K for the PA, TI, and CH. Inspection of the results at T = 100 K shows that in using PA, the partition function converges to within 1% of the exact partition function value of 1.69 106 at a k-level of 8 (results not shown). Using TI reduces this to k = 6, whereas with CH the desired accuracy is reached with k = 4 with the optimal t0.

Figure 2. Partition functions for the MO calculated by the CH algorithm at T = 100, 300, and 500 K, with varying values of the parameter t0.

Figure 3. Partition functions for the SDW potential calculated by the CH algorithm at T = 100, 300, and 500 K, with varying values of the parameter t0.

Table 3. Partition Functions for the MO Calculated by the PA, TI, and CH Algorithms at Various k-Levels at Temperatures T = 100, 300, and 500 K Q (T = 100 K) log2 P

PA

1

1.177 103

1.639 105

2.200 10

5

6.363 10

8

1.133 10

7

1.501 10

10

2 3

a

Q (T = 300 K) a

TI

CH

PA

1.491 109 12

6.639 10

13

5.354 10

Q (T = 500 K) a

TI

CH

PA

CHa

TI

1.046 102

1.201 103

1.019 104

2.818 102

7.739 103

3.317 103

3

2.046 10

4

5

3

3.963 10

3

3.140 103

8.615 10

5

3.240 10

3

3.126 103

1.574 10

4

2.938 10

7.054 10

5

6.704 10

9.571 10

3

4.741 10

4 5

10

3.978 10 5.728 1012

12

2.068 10 4.235 1013

13

3.156 10 2.989 1013

4

1.115 10 7.721 105

5

6.906 10 6.698 105

5

6.681 10 6.679 105

3

3.533 10 3.228 103

3

3.135 10 3.126 103

3.125 103 3.125 103

6

7.519 1013

3.093 1013

2.980 1013

6.935 105

6.681 105

6.679 105

3.151 103

3.125 103

3.125 103

The CH algorithm employed t0 = 0.1266. 1280



ARTICLE

Table 4. Partition Functions for the SDW Potential Calculated by the PA, TI, and CH Algorithms at Various k-Levels at Temperatures T = 100, 300, and 500 K Q (T = 100 K)

Q (T = 500 K) a

CHa

log2 P

PA

1

3.490 102

4.542 103

2.885 103

2.524 101

1.889 101

1.754 101

5.225 101

4.693 101

4.589 101

2

3

3

1

1

1

1

1

4.587 101

2

TI

1.103 10

CH

3.122 10

PA

2.649 10

TI

2.015 10

CH

1.780 10

PA

1.747 10

TI

4.788 10

4.601 10

3 4

3

4.879 10 3.208 103

3

2.838 10 2.618 103

3

2.577 10 2.570 103

1

1.826 10 1.768 101

1

1.751 10 1.747 101

1

1.747 10 1.747 101

1

4.642 10 4.601 101

1

4.588 10 4.587 101

4.587 101 4.587 101

5

2.740 103

2.574 103

2.570 103

1.752 101

1.747 101

1.747 101

4.591 101

4.587 101

4.587 101

3

3

3

1

1

1

1

1

4.587 101

6 a

Q (T = 300 K) a

2.613 10

2.570 10

2.570 10

1.748 10

1.747 10

1.747 10

4.588 10

4.587 10

The CH algorithm employed t0 = 0.1266.

Figure 4. Partition functions for the ADW potential calculated by the CH algorithm at T = 100, 300, and 500 K, with varying values of the parameter t0.

At T = 300 K the exact partition function value is 1.19 102. PA converges with k = 6, TI requires k = 4, and CH converges with k = 2. At T = 500 K the exact value of the partition function is 7.10 102. To reach convergence PA requires k = 5, TI requires k = 3, whereas CH needs k = 1. Thus, at all temperatures, the CH method yields a 25% enhancement in performance compared to TI. Moreover, CH reaches the performance of PA at ca. 20% of the cost at all temperatures. 4.5. H2O Molecule. To investigate the application of the various potentials described in this paper, we employed a water molecule at the B3LYP/6-31þG(d,p) level of theory. Specifically, we compute the quantum correction at T = 300 K, where the quantum effects are modest. These simulations show that the

CH potential is equally applicable to more complex potentials (Figure 5 and Table 6). At this temperature, TI and CH perform similarly well, while PA requires approximately four times as many beads. 4.6. Isotope Effect on KetoEnol Tautomerism in Alanine Racemase (AlaR). PLP is an essential cofactor for ubiquitous enzyme catalyzed transformations of amines and amino acids, such as racemizations, transaminations, and decarboxylations. A crucial question in all PLP-dependent enzymes is the tautomeric nature of the Schiff-base (Scheme 1), as it may exist in either the iminophenoxide or the enolimine form. The tautomeric form and hence the Schiff-base hydrogen-bond strength is highly sensitive to solvent polarity, and this topic has been addressed experimentally by NMR studies of the hydrogen-bond EIE.7577 From a computational perspective, it is therefore important to develop methods which can accurately predict the hydrogenbond EIE in enzymes. Herein, we employ the PA, TI, and CH methods with the mass-perturbation staging algorithm derived in eqs 3138. Initially, to validate the vibrational frequencies of the OH and NH stretches in tautomers of the pyridoxal moiety, we performed model calculations on the tautomeric ketoamine NH and enolimine OH forms (Scheme 1). The computed vibrational frequency for the NH-stretch in zwitterionic ketoamine NH tautomer was 3211.7 and 3290.9 cm1 at the target M06/631þG(d,p) level78 and at the AM1-SRP level, respectively (Table 7). The computed vibrational frequency for the OHstretch in the nonzwitterionic enolimine OH tautomer was 3413.6 and 3435.9 cm1 at the target M06/6-31þG(d,p) level and at the AM1-SRP level, respectively. Thus, the differences between the vibrational frequencies of the two tautomeric forms are 201.9 and 145.1 cm1 at the M06 and semiempirical levels. The computed gas-phase equilibrium IE is 0.89 and 1.06 at the M06 and semiempirical levels, respectively, where we have employed a scaling factor of 0.98 for M06 frequencies.79 The enzyme quantum simulations employed the PA, TI, and CH approaches in conjunction with the mass-perturbation staging algorithm. The PLP cofactor in AlaR is presented in Figure 6 with a protonated and deuterated Schiff base, and the numerical results obtained using eqs 1 and 38 are summarized in Tables 8 and 9. We estimate the converged value of the EIE at 298 K as 1.16 ( 0.06. Interestingly, PA shows robust performance with the mass-perturbation staging algorithm for computation of the EIE. Using only 3 or 6 beads, the EIE is estimated to be 1.19. With 12 beads, the result is 1.16, while further increasing the number of beads to 24 or 48 yields 1.15. Surprisingly, the TI and CH show poor performance using 3 or 6 beads. Using TI the 1281



ARTICLE

Table 5. Partition functions for the ADW potential calculated by the PA, TI, and CH algorithms at various k-levels at temperatures T = 100, 300, and 500 K Q (T = 100 K)

Q (T = 500 K) a

CHa

log2 P

PA

1

5.246 103

3.091 104

4.918 106

4.501 102

1.775 102

1.215 102

1.125 101

7.711 102

7.108 102

4

5

6

2

2

2

2

2

7.097 102

2

TI

4.274 10

CH

1.858 10

PA

2.065 10

TI

2.162 10

CH

1.303 10

PA

1.193 10

TI

8.375 10

7.189 10

3 4

5

3.273 10 5.462 106

6

3.312 10 1.916 106

6

1.725 10 1.693 106

2

1.457 10 1.261 102

2

1.206 10 1.193 102

2

1.192 10 1.192 102

2

7.446 10 7.186 102

2

7.105 10 7.097 102

7.096 102 7.096 102

5

2.427 106

1.714 106

1.691 106

1.209 102

1.192 102

1.192 102

7.119 102

7.096 102

7.096 102

6

6

6

2

2

2

2

2

7.096 102

6 a

Q (T = 300 K) a

1.863 10

1.692 10

1.691 10

1.196 10

1.192 10

1.192 10

7.102 10

7.096 10

The CH algorithm employed t0 = 0.1266.

Figure 5. Quantized water molecule treated at the B3LYP/6-31þG(d, p) level with 18 beads (P = 18).

Table 6. Vibrational Quantum f Classical Free Energy Correctiona (kcal/mol) for a Water Molecule Treated at the B3LYP/6-31þG(d,p) Level Computed Using the Staging Algorithm in Conjunction with the PA, TI, and CH Algorithms with Various Numbers of Beads (P) at Temperature T = 300 K P

PA

TI

CH

3

4.09 ( 0.06

7.38 ( 0.94

7.25 ( 0.48

6 12

6.70 ( 0.14 8.42 ( 0.22

8.89 ( 0.65 9.39 ( 0.31

8.86 ( 0.46 9.40 ( 0.22

24

9.08 ( 0.17

9.40 ( 0.26

9.27 ( 0.20

a

The calculations used eq 30 with 100 MC staging algorithm steps per classical point. The total number of classical points was 100. All values are averaged over 10 independent runs.

Table 7. Computed Unscaled Vibrational Frequencies (cm1) of the Schiff-Base Moiety in Tautomers of a Model Pyridoxal Compound M06/6-31þG(d,p)

AM1-SRP

ketoamine NH-form

3211.7

3290.9

enolimine OH-form

3413. 6

3435.9

EIE is estimated to be 1.58 and 1.25 using 3 and 6 beads, respectively. Using CH the EIE is estimated to be 1.44 and 1.19 using 3 and 6 beads, respectively. Further increasing the number of beads to 12, 24, or 48 yields converged values for both TI and

Figure 6. Quantized Schiff base and oxyanion in the PLP cofactor in AlaR with 48 beads and (A) protonated (B) deuterated Schiff base. The bead distributions were computed with the mass-perturbation staging algorithm.

CH. At the current simulation temperature and all numbers of beads employed, TI and CH are still expected to perform better than PA (i.e., OTI/CH(τ4) < OPA(τ2)). The reason for this seems 1282



ARTICLE

Table 8. RS and PS IE Computed for the KetoEnol Tautomerism in AlaR at T = 298 Ka RS P

PA

PS

TI

CH

PA

TI

CH

3

0.227 ( 0.009

0.039 ( 0.005

0.041 ( 0.005

0.270 ( 0.009

0.062 ( 0.005

0.059 ( 0.004

6

0.096 ( 0.004

0.034 ( 0.003

0.035 ( 0.003

0.114 ( 0.004

0.043 ( 0.002

0.041 ( 0.003

12

0.060 ( 0.003

0.042 ( 0.002

0.041 ( 0.002

0.070 ( 0.003

0.048 ( 0.002

0.047 ( 0.002

24

0.051 ( 0.002

0.046 ( 0.002

0.045 ( 0.002

0.059 ( 0.002

0.054 ( 0.002

0.053 ( 0.002

48

0.049 ( 0.002

0.047 ( 0.002

0.047 ( 0.002

0.056 ( 0.002

0.055 ( 0.002

0.054 ( 0.002

a

The calculations used eq 38 with 10 MC staging algorithm steps per classical point per isotope. The total number of classical points was 10 400. The 2 1/2 /N), which is the standard deviation in computing eq 38 in the RS or PS wells using N discrete points, and σi is the error is estimated as hσ = (∑N i=1σi ) standard deviation in computing IE at a discrete point in the reactant or product well.

Table 9. EIE Computed for the KetoEnol Tautomerism in AlaR at T = 298 Ka P

PA

TI

CH

3

1.19 ( 0.06

1.58 ( 0.25

1.44 ( 0.22

6

1.19 ( 0.07

1.25 ( 0.13

1.19 ( 0.13

12

1.16 ( 0.07

1.14 ( 0.08

1.15 ( 0.08

24 48

1.15 ( 0.06 1.15 ( 0.06

1.17 ( 0.07 1.15 ( 0.07

1.16 ( 0.07 1.16 ( 0.06

a

The calculations used eqs 1 and 38 with 10 MC staging algorithm steps per classical point per isotope. The total number of classical points was 2 2 1/2 10 400. The error is estimated as σ = ((σ hRS/IERS) þ (σ hPS/IEPS) ) 3 (IEPS/IERS), where IE = QL/QH is the IE in either the RS or PS wells, N 2 1/2 σ h = (∑i=1σi ) /N is the standard deviation in computing eq 38 in the RS or PS wells using N discrete points, and σi is the standard deviation in computing IE at a discrete point in the reactant or product well.

to be the greater uncertainty in the computed EIE using TI and CH with a low number of beads, which is interestingly due to the enhanced accuracy of the methods. Inspection of Table 8 reveals that the simulation error is reduced as the number of beads is increased. Moreover, the standard deviation is similar for all three methods. It is important to note that the standard deviation is not due to the sampling of the kinetic energy term, as this part is sampled exactly by the free-particle mass-perturbation staging algorithm, but rather due to sampling of the potential energy surface. Specifically, both the classical averaging over the potential energy surface as well as the PI sampling of the potential surface contribute to the standard deviation. This is due to the fluctuating nature of the complex potential energy surface in enzymes. Using TI and CH in computing eq 38 in the RS and PS, respectively, with a small number of beads yields fairly converged IE values with respect to number of beads (Table 8). In computing the EIE we need to divide IE in the PS and RS (EIE = IEPS/IERS), which in the case of TI and CH are small numbers with large error bars, yielding greater errors in the computed EIE (Table 9). On the other hand when using PA, the absolute error in computing the IE is greater due to the small number of beads and lack of higher order terms as in TI and CH. Thus, although PA exhibits greater absolute errors in computing the IE, these errors largely cancel out in the RS and PS, as the error is not in the leading digits of the IE. Finally, we compare the efficiency of the mass perturbation treatment using the staging and bisection algorithms. Specifically we computed the EIE in AlaR using 8 beads with the two methods, using either 10 or 20 MC steps per classical

configuration for a total of 5200 classical configurations. Employing 10 MC steps we obtained 1.19 ( 0.07 and 1.23 ( 0.07 for the staging and bisection algorithms, respectively, while using 20 MC steps we obtained 1.20 ( 0.05 for both the staging and bisection algorithms, respectively. Thus, the two methods give comparable results, and this conclusion is not expected to change when using a greater number of beads. Indeed, both the bisection and staging algorithms sample the kinetic part of the action exactly, and therefore for free particle sampling, their performance will be comparable. Thus, within the framework of eq 30 both sampling schemes may in principle be employed. However, the Chin action requires that the number of beads be a multiple of three (see eq 22) and may be readily achieved with the staging algorithm, which can sample any number of beads. However, the bisection algorithm naturally samples 2k number of beads in a naïve implementation, where k is the sampling level, and therefore is not generally suitable for the Chin action. Thus, the staging algorithm may be more flexible with respect to number of beads. We note that when the sampling entails not only the kinetic part of the action but also the potential part, the bisection algorithm may be advantageous. Using the bisection algorithm when moving P beads, the largest bead move is performed at the first MC step (i.e., the middle bead), and one may reject the collective move of P beads based on the move of a single bead (i. e., a single energy and force calculation as opposed to P such calculations).

5. DISCUSSION In this study we initially compare the performance of the PA algorithm and the higher order TI and CH algorithms on four model potentials: HO, MO, SDW, and ADW. We find that the CH algorithm with optimal parameters performs considerably better than the PA when computing the partition function for the model chemical potentials. This conclusion is in accordance with the findings of Sakkos et al for quantum liquids.58 The use of the CH approach is most beneficial at low temperatures where quantum effects are more pronounced. Nonetheless, we find that the TI approach performs nearly as well as CH, and the main gain is in going beyond the PA. These findings for model systems are of great importance when moving to condensed phase systems, where the addition of numerical noise complicates the performance analysis of the methods. The parametrized CH algorithm is expected to be of value in condensed phase simulations where the computational bottleneck is the energy evaluation, such as in simulations employing fully QM or hybrid QM/MM potential energy surfaces. In typical 1283


Journal of Chemical Theory and Computation uses of such potentials, iterative self-consistent field calculations are required in evaluating the potential energy, and these are computationally expensive. The computation of gradients, on the other hand, requires less effort than the energy evaluation itself. Thus, a PI method which can significantly reduce the number of energy evaluations is of great value. Indeed, the efficiency of the CH factorization in the calculation of nuclear QM effects using a complicated potential energy surface is exemplified in this work by calculations of water treated with DFT. This conclusion is also correct for a considerably more complex hybrid QM/MM potential energy surface such as the one employed here in the case of the enzyme AlaR. However, in computing IEs, numerical noise hampers the performance of both TI and CH with a small number of beads although the quantum effects are treated more accurately than with PA. This is largely due to the simulation noise inherent to any sampling method and not due to inherent properties of TI or CH. This is clear from the model calculations of the IE for the HO, where TI and CH displayed superb performance. Remarkably, PA is highly accurate in computing the EIE on the tautomerism of PLP in AlaR, even when using only three beads, when employing the mass-perturbation staging algorithm. The approaches employed in this work (eqs 30 and 38) are equally applicable to computing the centroid potential of mean force, and this is currently being investigated in our group. Additionally, PI schemes based on the flux autocorrelation methods which require the calculation of the entire density matrix will benefit from the CH algorithm. In such chemical rate calculations, the PA and TI approaches are expected to be much less efficient. Higher dimensionality derivations of flux autocorrelation methods in conjunction with the CH method are being pursued in our group. It is interesting to note that the enzyme environment enhances the EIE when compared to the gas-phase results. In the gas-phase, the EIE is computed to be 1.06, while in AlaR it is estimated as 1.16. This is indicative of a weakening of the intramolecular hydrogen bond relative to the gas phase. This is indeed expected as the highly polar active site in AlaR reduces the difference between the zero-point energies of the iminophenoxide or enolimine forms. This is in agreement with experimental work on model PLP systems in solvents of varying degrees of polarity.77 We believe the current approach will be of great use in the study of the effect of active site polarity on the hydrogen-bond strength in PLPdependent enzymes as well as other enzymes. Finally, it may be instructive to compare the current approach for computing IEs to other related approaches. Recently, Wong et al. employed classical TST, PI quantum TST, and the quantum instanton approaches to evaluate the quantized potential of mean force and KIE in malonaldehyde.80 In this study, the latter two approaches were found to give KIEs in reasonable agreement with each other, although a clear relationship between the two methods has not yet been established. The current mass-perturbation staging approach is in principle similar to the PI quantum TST with thermodynamic integration approach employed by Wong et al. However, the advantage of the current approach is that PI sampling is only required for the light and heavy isotopes, and no sampling of intermediate mass values is required during the perturbation. This one-step perturbation is highly efficient and is possible due to the fact that the current PI sampling is performed using free-particle MC.

6. CONCLUSIONS Higher-order corrections to the primitive approximation (PA) may considerably enhance the performance of quantum

ARTICLE

simulation methods. In this report we compare the PA and the higher order TakahashiImada (TI) algorithm with the gradient-based forward corrector algorithm due to Chin (CH) on a variety of model potentials. We find a unique parameter for the Chin algorithm which gives a good performance for all model potentials tested. Moreover, the PA, TI, and CH factorizations are employed to compute the quantum correction to a water molecule treated with the B3LYP functional with a 6-31þG(d,p) basis set. Finally, we employ the PA, TI, and CH methods to compute the equilibrium IE on the Schiff baseoxyanion tautomerism in the cofactor pyridoxal-50 -phosphate in the enzyme alanine racemase using a novel mass-perturbation staging algorithm. We find that the Chin algorithm performs well for the complex molecular systems as well, although numerical noise might hamper its performance when computing IEs with a small number of beads.

’ ASSOCIATED CONTENT

bS

Supporting Information. Figures of model potentials at temperatures 100500 K. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work has been supported by a start-up grant from BarIlan University, Israel Science Foundation, and the Alon fellowship from the Council for higher education, planning, and budgeting committee. We thank Prof. Siu A. Chin for helpful comments to this manuscript. ’ REFERENCES (1) Frey, P. A.; Hegeman, A. D. Enzymatic reaction mechanisms; Oxford University Press: New York, 2007. (2) Kohen, A.; Limbach, H. H. Isotope Effects in Chemistry and Biology; Taylor and Francis Group, CRC Press: New York, 2006. (3) Lewis, B. E.; Schramm, V. L. Enzymatic binding isotope effects and the interaction of glucose with hexokinase. In Isotope Effects in Chemistry and Biology; Kohen, A., Limbach, H. H., Eds.; Taylor and Francis Group, CRC Press: New York, 2006; pp 1019. (4) Limbach, H. H.; Denisov, G. S.; Golubev, N. S., Hydrogen bond isotope effects studied by NMR. In Isotope Effects in Chemistry and Biology; Kohen, A., Limbach, H. H., Eds.; Taylor and Francis Group, CRC Press: New York, 2006; pp 193. (5) Voth, G. A. Feynman path integral formulation of quantum mechanical transition state theory. J. Phys. Chem. 1993, 97, 8365. (6) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. (7) Berne, B. J.; Thirumalai, D. On the simulation of quantumsystems - path integral methods. Annu. Rev. Phys. Chem. 1986, 37, 401. (8) Marx, D.; Parrinello, M. Structural quantum effects and 3-center 2-electron bonding in ch5þ. Nature 1995, 375, 216. (9) Hwang, J. K.; Chu, Z. T.; Yadav, A.; Warshel, A. Simulations of quantum-mechanical corrections for rate constants of hydride-transfer reactions in enzymes and solutions. J. Phys. Chem. 1991, 95, 8445. (10) Hwang, J. K.; Warshel, A. A quantized classical path approach for calculations of quantum-mechanical rate constants. J. Phys. Chem. 1993, 97, 10053. 1284


Journal of Chemical Theory and Computation (11) Major, D. T.; Gao, J. L. Implementation of the bisection sampling method in path integral simulations. J. Mol. Graphics Modell. 2005, 24, 121. (12) Major, D. T.; Garcia-Viloca, M.; Gao, J. L. Path integral simulations of proton transfer reactions in aqueous solution using combined QM/MM potentials. J. Chem. Theory Comput. 2006, 2, 236. (13) Major, D. T.; Gao, J. L. An integrated path integral and freeenergy perturbation-umbrella sampling method for computing kinetic isotope effects of chemical reactions in solution and in enzymes. J. Chem. Theory Comput. 2007, 3, 949. (14) Gao, J. L.; Wong, K. Y.; Major, D. T. Combined QM/MM and path integral simulations of kinetic isotope effects in the proton transfer reaction between nitroethane and acetate ion in water. J. Comput. Chem. 2008, 29, 514. (15) Wang, M. L.; Lu, Z. Y.; Yang, W. T. Nuclear quantum effects on an enzyme-catalyzed reaction with reaction path potential: Proton transfer in triosephosphate isomerase. J. Chem. Phys. 2006, 124, 124516. (16) Wang, Q.; Hammes-Schiffer, S. Hybrid quantum/classical path integral approach for simulation of hydrogen transfer reactions in enzymes. J. Chem. Phys. 2006, 125, 184102. (17) Marx, D.; Tuckerman, M. E.; Martyna, G. J. Quantum dynamics via adiabatic ab initio centroid molecular dynamics. Comput. Phys. Commun. 1999, 118, 166. (18) Tuckerman, M. E.; Marx, D.; Klein, M. L.; Parrinello, M. On the quantum nature of the shared proton in hydrogen bonds. Science 1997, 275, 817. (19) Tuckerman, M. E.; Marx, D.; Parrinello, M. The nature and transport mechanism of hydrated hydroxide ions in aqueous solution. Nature 2002, 417, 925. (20) Paesani, F.; Iuchi, S.; Voth, G. A. Quantum effects in liquid water from an ab initio-based polarizable force field. J. Chem. Phys. 2007, 127, 074506. (21) Ohta, Y.; Ohta, K.; Kinugawa, K. Quantum effect on the internal proton transfer and structural fluctuation in the H-5(þ) cluster. J. Chem. Phys. 2004, 121, 10991. (22) Hayashi, A.; Shiga, M.; Tachikawa, M. H/D isotope effect on the dihydrogen bond of NH4þ center dot center dot center dot BeH2 by ab initio path integral molecular dynamics simulation. J. Chem. Phys. 2006, 125, 204310. (23) Wong, K. Y.; Gao, J. An automated integration-free pathintegral method based on Kleinert’s variational perturbation theory. J. Chem. Phys. 2007, 127, 211103. (24) Wong, K. Y.; Gao, J. Systematic approach for computing zeropoint energy, quantum partition function, and tunneling effect based on Kleinert’s variational perturbation. J. Chem. Theory Comput. 2008, 4, 1409. (25) Alhambra, C.; Corchado, J.; Sanchez, M. L.; Garcia-Viloca, M.; Gao, J.; Truhlar, D. G. Canonical variational theory for enzyme kinetics with the protein mean force and multidimensional quantum mechanical tunneling dynamics. Theory and application to liver alcohol dehydrogenase. J. Phys. Chem. B 2001, 105, 11326. (26) Pu, J. Z.; Gao, J. L.; Truhlar, D. G. Multidimensional tunneling, recrossing, and the transmission coefficient for enzymatic reactions. Chem. Rev. 2006, 106, 3140. (27) Billeter, S. R.; Webb, S. P.; Agarwal, P. K.; Iordanov, T.; Hammes-Schiffer, S. Hydride transfer in liver alcohol dehydrogenase: Quantum dynamics, kinetic isotope effects, and role of enzyme motion. J. Am. Chem. Soc. 2001, 123, 11262. (28) Iyengar, S. S.; Sumner, I.; Jakowski, J. Hydrogen tunneling in an enzyme active site: A quantum wavepacket dynamical perspective. J. Phys. Chem. B 2008, 112, 7601. (29) Antoniou, D.; Basner, J.; Nunez, S.; Schwartz, S. D. Computational and theoretical methods to explore the relation between enzyme dynamics and catalysis. Chem. Rev. 2006, 106, 3170. (30) Major, D. T.; York, D. M.; Gao, J. L. Solvent polarization and kinetic isotope effects in nitroethane deprotonation and implications to the nitroalkane oxidase reaction. J. Am. Chem. Soc. 2005, 127, 16374.

ARTICLE

(31) Major, D. T.; Nam, K.; Gao, J. L. Transition state stabilization and alpha-amino carbon acidity in alanine racemase. J. Am. Chem. Soc. 2006, 128, 8114. (32) Major, D. T.; Gao, J. L. A combined quantum mechanical and molecular mechanical study of the reaction mechanism and alpha-amino acidity in alanine racemase. J. Am. Chem. Soc. 2006, 128, 16345. (33) Rubinstein, A.; Major, D. T. Catalyzing Racemizations in the Absence of a Cofactor: The Reaction Mechanism in Proline Racemase. J. Am. Chem. Soc. 2009, 131, 8513. (34) Major, D. T.; Heroux, A.; Orville, A. M.; Valley, M. P.; Fitzpatrick, P. F.; Gao, J. Differential quantum mechanical tunneling in the uncatalyzed and in the Nitroalkane Oxidase proton abstraction of nitroethane. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 20734. (35) Zimmermann, T.; Vanicek, J. Path integral evaluation of equilibrium isotope effects. J. Chem. Phys. 2009, 131, 024111. (36) Zimmermann, T.; Vanicek, J. Three applications of path integrals: equilibrium and kinetic isotope effects, and the temperature dependence of the rate constant of the [1,5] sigmatropic hydrogen shift in (Z)-1,3-pentadiene. J. Mol. Model. 2010, 16, 1779. (37) Sprik, M.; Klein, M. L.; Chandler, D. Staging - a sampling technique for the Monte Carlo evaluation of path-integrals. Phys. Rev. B: Condens. Matter Mater. Phys. 1985, 31, 4234. (38) Pollock, E. L.; Ceperley, D. M. Simulation of quantum manybody systems by path-integral methods. Phys. Rev. B: Condens. Matter Mater. Phys. 1984, 30, 2555. (39) Cao, J.; Berne, B. J. A Born-Oppenheimer approximation for path integrals with an application to electron solvation in polarizable fluids. J. Chem. Phys. 1993, 99, 2902. (40) Takahashi, M.; Imada, M. Monte Carlo calculation of quantumsystems. 2. Higher-order correction. J. Phys. Soc. Jpn. 1984, 53, 3765. (41) Li, X. P.; Broughton, J. Q. High-order correction to the Trotter expansion for use in computer simulation. J. Chem. Phys. 1987, 86, 5094. (42) Suzuki, M. Compact exponential product formulas and operator functional derivative. J. Math. Phys. 1997, 38, 1183. (43) Suzuki, M. Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations. Phys. Lett. A 1990, 146, 319. (44) Chin, S. A. Symplectic integrators from composite operator factorizations. Phys. Lett. A 1997, 226, 344. (45) Chin, S. A.; Chen, C. R. Fourth order gradient symplectic integrator methods for solving the time-dependent Schrodinger equation. J. Chem. Phys. 2001, 114, 7338. (46) Chin, S. A.; Chen, C. R. Gradient symplectic algorithms for solving the Schrodinger equation with time-dependent potentials. J. Chem. Phys. 2002, 117, 1409. (47) Forbert, H. A.; Chin, S. A. Fourth-order diffusion Monte Carlo algorithms for solving quantum many-body problems. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 63, 144518. (48) De Raedt, H.; De Raedt, B. Applications of the generalized Trotter formula. Phys. Rev. A: At., Mol., Opt. Phys. 1983, 28, 3575. (49) Kono, H.; Takasaka, A.; Lin, S. H. Monte Carlo calculation of the quantum partition function via path integral formulations. J. Chem. Phys. 1988, 88, 6390. (50) Schwartz, S. D. Accurate quantum mechanics from high order resummed operator expansions. J. Chem. Phys. 1994, 100, 8795. (51) Weht, R. O.; Kohanoff, J.; Estrin, D. A.; Chakravarty, C. An ab initio path integral Monte Carlo simulation method for molecules and clusters: Application to Li-4 and Li-5(þ). J. Chem. Phys. 1998, 108, 8848. (52) Jang, S. J.; Jang, S. M.; Voth, G. A. Applications of higher order composite factorization schemes in imaginary time path integral simulations. J. Chem. Phys. 2001, 115, 7832. (53) Omelyan, I. P.; Mryglod, I. M.; Folk, R. Construction of highorder force-gradient algorithms for integration of motion in classical and quantum systems. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2002, 66, 026701. (54) Yamamoto, T. M. Path-integral virial estimator based on the scaling of fluctuation coordinates: Application to quantum clusters with fourth-order propagators. J. Chem. Phys. 2005, 123, 104101. 1285


Journal of Chemical Theory and Computation (55) Whitfield, T. W.; Martyna, G. J. Low variance energy estimators for systems of quantum Drude oscillators: Treating harmonic path integrals with large separations of time scales. J. Chem. Phys. 2007, 126, 074104. (56) Cuervo, J. E.; Roy, P. N.; Boninsegni, M. Path integral ground state with a fourth-order propagator: Application to condensed helium. J. Chem. Phys. 2005, 122, 114504. (57) Chin, S. A. Quantum statistical calculations and symplectic corrector algorithms. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2004, 69, 046118. (58) Sakkos, K.; Casulleras, J.; Boronat, J. High order Chin actions in path integral Monte Carlo. J. Chem. Phys. 2009, 130, 204109. (59) Klemm, A. D.; Storer, R. G. Structure of quantum fluids helium and neon. Aust. J. Phys. 1973, 26, 43. (60) Thirumalai, D.; Berne, B. J. On the calculation of time correlation-functions in quantum-systems - path integral techniques. J. Chem. Phys. 1983, 79, 5029. (61) Thirumalai, D.; Bruskin, E. J.; Berne, B. J. An iterative scheme for the evaluation of discretized path-integrals. J. Chem. Phys. 1983, 79, 5063. (62) Sethia, A.; Sanyal, S.; Singh, Y. Discretized path integral method and properties of a quantum system. J. Chem. Phys. 1990, 93, 7268. (63) Chandler, D.; Wolynes, P. G. Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids. J. Chem. Phys. 1981, 74, 4078. (64) Trotter, H. F. On the product of semi-groups of operators. Proc. Am. Math. Soc. 1959, 10, 545. (65) Scuro, S. R.; Chin, S. A. Forward symplectic integrators and the long-time phase error in periodic motions. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2005, 71, 056703. (66) Khristenko, S. V.; Maslov, A. I.; Shevelko, V. P. Molecules and Their Spectroscopic Properties; Springer: Berlin, Heidelberg, Germany, 1998. (67) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A Program for Macromolecular Energy, Minimization, and Dynamics Calculations. J. Comput. Chem. 1983, 4, 187. (68) Becke, A. D. Density-functional thermochemistry III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648. (69) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785. (70) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; John Wiley & Sons: New York, 1986. (71) Guest, M. F.; Bush, I. J.; van Dam, H. J. J.; Sherwood, P.; Thomas, J. M. H.; van Lenthe, J. H.; Havenith, R. W. A.; Kendrick, J. The GAMESS-UK electronic structure package: algorithms, developments and applications. Mol. Phys. 2005, 103, 719. (72) Mielke, S. L.; Truhlar, D. G. Displaced-points path integral method for including quantum effects in the Monte Carlo evaluation of free energies. J. Chem. Phys. 2001, 115, 652. (73) Morse, P. M. Diatomic molecules according to the wave mechanics. II. Vibrational levels. Phys. Rev. 1929, 34, 57. (74) Morrone, J. A.; Lin, L.; Car, R. Tunneling and delocalization effects in hydrogen bonded systems: A study in position and momentum space. J. Chem. Phys. 2009, 130, 204511. (75) Sharif, S.; Denisov, G. S.; Toney, M. D.; Limbach, H. H. NMR studies of solvent-assisted proton transfer in a biologically relevant Schiff base: Toward a distinction of geometric and equilibrium H-bond isotope effects. J. Am. Chem. Soc. 2006, 128, 3375. (76) Sharif, S.; Schagen, D.; Toney, M. D.; Limbach, H. H. Coupling of functional hydrogen bonds in pyridoxal-50 -phosphate-enzyme model systems observed by solid-state NMR spectroscopy. J. Am. Chem. Soc. 2007, 129, 4440. (77) Sharif, S.; Denisov, G. S.; Toney, M. D.; Limbach, H. H. NMR studies of coupled low- and high-barrier hydrogen bonds in pyridoxal-50 phosphate model systems in polar solution. J. Am. Chem. Soc. 2007, 129, 6313.

ARTICLE

(78) Zhao, Y.; Truhlar, D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215. (79) Alecu, I. M.; Zheng, J.; Zhao, Y.; Truhlar, D. G. Computational thermochemistry: Scale factor databases and scale factors for vibrational frequencies obtained from electronic model chemistries. J. Chem. Theory Comput. 2010, 6, 2872. (80) Wong, K. F.; Sonnenberg, J. L.; Paesani, F.; Yamamoto, T.; Vanicek, J.; Zhang, W.; Schlegel, H. B.; Case, D. A.; Cheatham, T. E., III; Miller, W. H.; Voth, G. A. Proton transfer studied using a combined ab initio reactive potential energy surface with quantum path integral methodology. J. Chem. Theory. Comp. 2010, 6, 2566.

1286


Path-Integral Calculations of Nuclear Quantum Effects in Model

Recommend Documents