First-year progress report on NSF grant CHE03-49122

 

"Quantum Mechanical Effects in Complex Systems"

 

Principal investigator: Donald G. Truhlar

 

March 31, 2005

 

691. "Introductory Lecture: Nonadiabatic Effects in Chemical Dynamics," A. W. Jasper, C. Zhu, S. Nangia, and D. G. Truhlar, Faraday Discussions 127, 1-22 (2004).

 

Recent progress in the theoretical treatment of electronically nonadiabatic processes is discussed. First we discuss the generalized Born–Oppenheimer approximation, which identifies a subset of strongly coupled states. Then we discuss the relative advantages and disadvantages of adiabatic and diabatic representations of the coupled surfaces and their interactions. Ab initio diabatic representations that do not require tracking geometric phases or calculating singular nonadiabatic nuclear momentum coupling are presented as one promising approach for characterizing the coupled electronic states of polyatomic photochemical systems. Such representations can be accomplished by methods based on functionals of the adiabatic electronic density matrix and the identification of reference orbitals for use in an overlap criterion. Next, four approaches to calculating or modeling electronically nonadiabatic dynamics are discussed: (1) accurate quantum mechanical scattering calculations, (2) approximate wave packet methods, (3) surface hopping, and (4) self-consistent-potential semiclassical approaches. The last two of these are particularly useful for polyatomic photochemistry, and recent refinements of these approaches are discussed. For example, considerable progress has been achieved in making the surface hopping method more applicable to the study of systems with weakly coupled electronic states. This includes introducing uncertainty principle considerations to alleviate the problem of classically forbidden surface hops and the development of an efficient sampling algorithm for low-probability events. A topic whose central importance in a number of quantum mechanical fields is becoming more widely appreciated is the introduction of decoherence into the quantal degrees of freedom to account for the effect of the classical treatment on the other degrees of freedom, and we discuss how the introduction of such decoherence into a self-consistent-potential approximation leads to a reasonably accurate but very practical trajectory method for electronically nonadiabatic processes. Finally, the performances of several dynamical methods for Landau–Zener-type and Rosen–Zener–Demkov-type reactive scattering problems are compared.

 

687. "Coherent Switching with Decay of Mixing: An Improved Treatment of Electronic Coherence for Non-Born-Oppenheimer Trajectories," C. Zhu, S. Nangia, A. W. Jasper, and D. G. Truhlar, Journal of Chemical Physics 121, 7658-7670 (2004).

 

The self-consistent decay-of-mixing (SCDM) semiclassical trajectory method for electronically nonadiabatic dynamics, developed under our previous NSF grant, was improved by modifying the switching probability that determines the instantaneous electronic state toward which the system decoheres, which is the pointer state in the language of quantum information theory. The method is called coherent switching with decay of mixing, CSDM, and it differs from the previously presented SCDM method in that the electronic amplitudes controlling the switching of the pointer state are treated fully coherently in the electronic equations of motion for each complete passage through a strong interaction region. It was tested against accurate quantum mechanical calculations for 12 atom-diatom scattering test cases. Also tested were the SCDM method and the trajectory surface hopping method of Parlant and Gislason that requires coherent passages through each strong interaction region, and which we call the ‘‘exact complete passage’’ trajectory surface hopping (ECP-TSH) method. The results are compared with previously presented results for the fewest switches with time uncertainty (FSTU, our own surface hopping method, developed under our previus NSF grant) and Tully’s fewest switches (TFS) surface hopping methods and the semiclassical Ehrenfest method. We find that the CSDM method is the most accurate of the semiclassical trajectory methods tested. Including coherent passages improves the accuracy of the SCDM method, such that the CSDM method is now the most accurate available method for quantum photochemistry simulations. Interestingly, though, the ECP-TSH method is not more accurate on average than the FSTU method, which remains the most accurate surface hopping method. We can understand this qualitative difference in terms of the more physical way in which coherent passages and reduction of the wave packet (reinitialization of the density matrix) are built into the CSDM method.

 

690. "Ensemble-Averaged Variational Transition State Theory with Optimized Multidimensional Tunneling for Enzyme Kinetics and Other Condensed-Phase Reactions," D. G. Truhlar, J. Gao, M. Garcia-Viloca, C. Alhambra, J. Corchado, M. L. Sanchez, and T. D. Poulsen, International Journal of Quantum Chemistry 100, 1136-1152 (2004). (Proceedings of the 2004 Sanibel Symposium)

 

Variational transition state theory was originally reduced to practice for reactions of small molecules in the gas phase (our gas-phase work on VTST is supported by the U.S. Department of Energy). In our previous NSF grant, we presented a new version of VTST, based on free energy sampling and semiclassical tunneling, that is suitable for complex reactions in condensed phases, in particular for the case where there is an ensemble of reaction paths and strong coupling of the reactive subsystem (primary subsystem) to its environment. A key advantage of our approach is the consistent inclusion of quantum mechanical effects, including validated methods developed in our group for mutltidimensional unneling. This paper provides an overview of the new method developed to include quantum mechanical effects and free energy sampling in calculations of reaction rates in enzymes. The paper includes an overview of variational transition state theory with optimized multidimensional tunneling for simple gas-phase reactions and then shows how this is extended to incorporate free energy effects and to include protein motions in the reaction coordinate by ensemble averaging. Finally we summarize recent very successful comparisons to experiment for primary and secondary kinetic isotope effects for proton and hydride transfer reactions catalyzed by enzymes, including one case where we predicted the secondary kinetic isotope effect before it was measured.

 

694. "Accurate Vibrational-Rotational Partition Functions and Standard-State Free Energy Values for H2O2 from Monte Carlo Path-Integral Calculations," V. A. Lynch, S. L. Mielke, and D. G. Truhlar, Journal of Chemical Physics 121, 5148-5162 (2004).

 

Accurate quantum mechanical partition functions and absolute free energies of hydrogen peroxide are determined using the realistic potential energy surface of Koput, Carter, and Handy for temperatures ranging from 300 to 2400 K by using Monte Carlo path integral calculations with new, efficient polyatomic importance sampling methods. The path centroids are sampled in Jacobi coordinates via a set of independent ziggurat schemes, which provides a new approach to importance sampling that may have widespread applicability. The calculations employed enhanced-same-path extrapolation of trapezoidal Trotter Fourier path integrals, and the paths were constructed using fast Fourier sine transforms. Importance sampling was also used in Fourier coefficient space, and adaptively optimized stratified sampling was added on top of the importance sampling in configuration space. The free energy values obtained from the path-integral calculations are compared to separable-mode approximations, to the Pitzer–Gwinn approximation, and to values in thermodynamic tables. Our calculations support the recently proposed revisions to the JANAF tables. This is the first example in the literature of a converged quantum mechanical partition coefficient for a molecule with a torsion, and it that respect it provides a milestone for the quantum mechanical treatment of anharmonic vibrations.

 

698. "Density-Functional and Hybrid-DFT SM5.43R Continuum Solvation Models for Aqueous and Organic Solvents," J. D. Thompson, C. J. Cramer, and D. G. Truhlar, Theoretical Chemistry Accounts 113, 107-131 (2005).

 

Hybrid density functional theory, which is a combined Hartree–Fock and density functional method, provides a simple but effective way to incorporate nonlocal exchange effects and the static and dynamical correlation energy into an orbital-based theory with affordable computational cost for many important problems of gas-phase chemistry. The inclusion of a reaction field representing an implicit solvent in a self-consistent hybrid density functional calculation provides an effective and efficient way to extend this approach to problems of liquid-phase chemistry. In previous work, supported by NSF, we have parameterized several models based on this approach, and in the present article, we present several new parameterizations based on implicit solvation models SM5.43 and SM5.43R. In particular, we extend the applicability of these solvation models to several combinations of the MPWX hybrid-density functional with various one-electron basis sets, where MPWX denotes a combination of Barone and Adamo’s modified version of Perdew and Wang’s exchange functional, Perdew and Wang’s correlation functional, and a percentage X of exact Hartree–Fock exchange. By carrying out parametrizations for variable X, we open the door to applications employing specific reaction parameters. SM5.43R parameter optimizations are presented for the MPWX/MIDI!, MPWX/MIDI!6D, and MPWX/ 6-31+G(d,p) combinations with X = 0 (i.e., pure density functional theory), 25, 42.8, and 60.6, and for MPWX/6-31G(d) and MPWX/6-31+G(d), with X = 0, 42.8, and 60.6; this constitutes a total of 18 new parameter sets. [Note that parameter optimizations using MPW25/6-31G(d) and MPW25/6-31+G(d) were carried out in a previous SM5.43R parameterization.] For each of the five basis sets, we found no significant loss in the accuracy of the model when parameters averaged over the four values of X are used instead of the parameters optimized for a specific value of X. Therefore for each of the five basis sets used here, the SM5.43R and SM5.43 models are defined to have a single parameter set that can be used for any value of X between 0 and 60.6. The new models yield accurate free energies of solvation for a broad range of solutes in both water and organic solvents. On the average, the mean-unsigned errors, as compared with those from experiment, of the free energies of solvation of neutral solutes range from 0.50 to 0.55 kcal/mol and those for ions range from 4.5 to 4.9 kcal/mol. Since the SM5.43R model computes the free energy of solvation as a sum of bulk-electrostatic and non-bulk-electrostatic contributions, it may be used for detailed analysis of the physical effects underlying a calculation of the free energy of solvation. Several calculations illustrating the partitioning of these contributions for a variety of solutes in n-hexadecane, 1-octanol, and water are presented.

 

699. "Conical Intersections and Semiclassical Trajectories: Comparison to Accurate Quantum Dynamics and Analyses of the Trajectories," A. W. Jasper and D. G. Truhlar, J. Chem. Phys. 122, 044101/1-16 (2005).

 

Our previous tests of semiclassical methods were carried out for systems withweakly coupled surfaces and avoided crossings. In this paper we extend these tests to systems with conical intersections. In particular, semiclassical trajectory methods are tested for electronically nonadiabatic systems with conical intersections. Five triatomic model systems are presented, and each system features two electronic states that intersect via a seam of conical intersections (CIs). Fully converged, full-dimensional (six dimensions) quantum mechanical scattering calculations are carried out (with our locally developed VP program based on variational methods we developed with earlier NF support)) for all five systems at energies that allow for electronic de-excitation via the seam of CIs; these provide a benchmark for ascertaining the validity of simpler methods that can be applied to larger systems. Several semiclassical trajectory methods are tested against the accurate quantum mechanical results. For four of the five model systems, the diabatic representation is the preferred (most accurate) representation for semiclassical trajectories, as correctly predicted by the Calaveras County criterion that we developed with earlier NSF support. Four surface hopping methods are tested and have overall relative errors of 40%–60%. The semiclassical Ehrenfest method has an overall error of 66%, and the self-consistent decay of mixing (SCDM) and coherent switches with decay of mixing (CSDM) methods are the most accurate methods overall with relative errors of ~32%. Furthermore, the CSDM method is less representation dependent than both the SCDM and the surface hopping methods, making it the preferred semiclassical trajectory method. Finally, taking advantage of the ability of semiclassical methods to produce a physical picture of the dynamics, the behavior of semiclassical trajectories near conical intersections is discussed.

 

704. "Quantum Mechanical Reaction Rate Constants by Vibrational Configuration Interaction. Application to the OH + H2 -> H2O + H Reaction as a Function of Temperature," A. Chakraborty and D. G. Truhlar, Proceedings of the National Academy of Sciences, U.S.A., accepted and available on-line at PNAS Early Edition. (Invited paper for Special Feature issue on Chemical Theory and Computation)

 

There has been a tendency in recent years to view time-dependent quantum mechanics as the only efficient way to calculate converged quantum dynamics for polyatomic systems. In this paper, we show the advantages of time-independent quantum mechanics. The thermal rate constant of the three-dimensional reaction of hydroxide radical with hydrogen molecule to produce water and a hydrogen atom was computed by using the flux autocorrelation function, with a time-independent square-integrable basis set. Two modes that actively participate in bond making and bond breaking were treated by using two-dimensional distributed Gaussian functions, and the remaining (nonreactive) modes were treated by using harmonic oscillator functions. The finite-basis eigenvalues and eigenvectors of the Hamiltonian were obtained by solving the resulting generalized eigenvalue equation, and the flux autocorrelation function for a dividing surface optimized in reduced-dimensionality calculations was represented in the basis formed by the eigenvectors of the Hamiltonian. The rate constant was obtained by integrating the flux autocorrelation function. The choice of the final time to which the integration is carried was determined by a plateau criterion. The potential energy surface was from Wu, Schatz, Lendvay, Fang, and Harding (WSLFH). The calculated thermalrate constants were compared with values calculated by wave packets, and excellent agreement was found; this is a rare example of two entirely different quantum dynamical methods agreeing for polyatomic dynamics. The success of the present calculations demonstrates that time-independent vibrational configuration interaction can be a very convenient way to calculate converged quantum mechanical rate constants, and it opens the possibility of calculating converged rate constants for much larger reactions than have been treated until now.

 

708. "Generalized Hybrid Orbital Method for Combined Quantum Mechanical and Molecular Mechanical Calculations Based on Density Functional Theory and Hybrid Density Functional Theory," J. Pu, J. Gao,and D. G. Truhlar, ChemPhysChem, in press. (invited paper for a Festschrift issue dedicated to M. Parrinello—it is a secret from Michele—don't tell him!)

 

The generalized hybrid orbital (GHO) method has previously been formulated for combining molecular mechanics with various levels of quantum mechanics, in particular semiempirical neglect-of diatomic differential overlap theory, ab initio Hartree-Fock theory, and self-consistent-charge density-functional tight binding theory. In our previous paper extending the GHO method to be compatible with the ab initio Hartree-Fock treatment of electronic structure, we showed how it is necessary to consider overlap effects much more carefully then in semiempirical calculations, but that this can indeed be done consistently. This also opened the door to using DFT and hybrid DFT in order to include electron correlation effects accurately, and in the present article we extend the GHO method to DFT and hybrid DFT with the generalized gradient approximation (denoted GHO-DFT and GHO-HDFT, respectively) using Gaussian-type orbitals as basis functions. In the proposed GHO-(H)DFT formalism, charge densities in auxiliary hybrid orbitals are included for calculating the total electron density. The orthonormality constraints involving the auxiliary Kohn-Sham orbitals are satisfied by carrying out the hybridization in terms of a set of Löwdin symmetrically orthogonalized atomic basis functions. Analytical gradients are formulated for GHO-(H)DFT by incorporating additional forces associated with GHO basis transformations. Scaling parameters are introduced for some of the one-electron integrals and are optimized to improve the geometry near the QM/MM boundary region. The GHO-(H)DFT method based on GGA methods (BLYP and mPWPW91) and HDFT methods (B3LYP, mPW1PW91, and MPW1K) is tested for geometries and atomic charges against a set of small molecules. The following quantities are tested: (1) the C-C stretch potential in ethane, (2) the torsional barrier for internal rotation around the central C-C bond in n-butane, (3) proton affinities for a set of alcohols, amines, thiols, and acids (4) the conformational energies of alanine dipeptide, and (5) the barrier height of the hydrogen atom transfer reaction between n-butane and n-butyl radical, where the reaction center is described at the MPW1K/6-31G(d) level of theory.

 

707. "Redistributed Charge and Dipole Schemes for Combined Quantum Mechanical and Molecular Mechanical Calculations," H. Lin and D. G. Truhlar, Journal of Physical Chemistry B, in press, accepted in January, 2005; proofs returned in March, 2005.

 

The combined quantum mechanical and molecular mechanical (QM/MM) method is a powerful tool for studying many chemical and biochemical processes such as enzyme reactions. A QM/MM model treats a relatively localized region (e.g., where bond breaking/forming or electronic excitation occur) with QM methods and includes the influence of the surroundings at the MM level. In some cases, such as the treatment of solvation, the boundary between the QM and MM subsystems is between solute and solvent molecules, and no covalent bond is cut. In many other cases, however, passing the boundary through covalent bonds is desirable, and special care is required to treat the boundary. In this study we demonstrated the importance of correctly handling the MM partial point charges at the QM/MM boundary, and in particular, we made progress in further developing QM/MM methods in two important aspects: (1) Two schemes, namely, the redistributed charge (RC) scheme and the redistributed charge and dipole (RCD) scheme, were introduced to handle link atoms in QM/MM calculations. In both schemes, the point charge at the MM boundary atom that is replaced by the link atom is redistributed to the midpoint of the bonds that connect the MM boundary atom and its neighboring MM atoms. These redistributed charges serve as classical mimics for the auxiliary orbitals associated with the MM host atom in the generalized hybrid orbital (GHO) method. In the RCD scheme, the dipoles of these bonds are preserved by further adjustment of the values of the redistributed charges. The treatments are justified as classical analogues of the QM description given by the GHO method. (2) The new methods were compared quantitatively to similar methods that were suggested by previous work, namely, a shifted-charge scheme and three eliminated charge schemes (including the default of Gaussian03). The comparisons were carried out for a series of molecules in terms of proton affinities and geometries. Point charges derived from various charge models were tested. The results demonstrate that it is critical to preserve charge and bond dipole and that it is important to use accurate MM point charges in QM/MM boundary treatments. The RCD scheme was further applied to study the H atom transfer reaction of methyl radical with 1-propanol. Various QM levels of theory were tested to demonstrate the generality of the methodology. Encouragingly, we found that the QM/MM calculations obtained a reaction energy, barrier height, saddle-point geometry, and imaginary frequency at the saddle point in quite good agreement with full QM calculations at the same level. Furthermore, analysis based on energy decomposition revealed the quantitatively similar interaction energies between the QM and the MM subsystems for the reactant, for the saddle point, and for the product. These interaction energies almost cancel each other energetically, resulting in negligibly small net effects on the reaction energy and barrier height. However, the charge distribution of the QM atoms is greatly affected by the polarization effect of the MM point charges. The QM/MM charge distribution agrees much better with full QM results than does the unpolarized charge distribution of the capped primary subsystem. The fact that the new RCD scheme works well is especially exciting in lieu of the fact that it is a completely general scheme that may be used with any quantum mechanical level (correlation method and basis set), with any molecular mechanics method that uses atom-centered point charges for the electrostatics, and with any electronic structure program that allows both positive and negative point charges. There are no new parameters, no pseudopotentials, and no integral scaling. The method is illustrated by applications with several different levels of QM theory (Hartree-Fock theory, second-order perturbation theory, hybrid density functional theory, and coupled cluster theory) and several choices for the molecular mechanics parameters.

 

710. "Small Temperature Dependence of the Kinetic Isotope Effect for the Hydride Transfer Reaction Catalyzed by Escherichia coli,Dihydrofolate Reductase," J. Pu, S. Ma, J. Gao, and D. G. Truhlar, Journal of Physical Chemistry B, accepted as a Letter on April 1, 2005. Manuscript no. JP051184C

 

Many enzyme reactions involve transfer of a hydron (proton, hydride ion, or hydrogen atom); such reactions have significant contributions from zero point energy and tunneling. Kinetic isotope effects have been extensively used to study these effects. Recent multidimensional tunneling calculations have been successfully used to study tunneling effects on enzymatic KIEs, but the most accurate methods have so far been applied only at a single temperature. Surprisingly, a number of enzymes were found to display almost temperature-independent KIEs, which are contrary to experience with small-molecule chemistry or simple tunneling models. This paper presents a theoretical explanation of the unusual temperature dependence of kinetic isotope effects (KIEs) for the hydride transfer reaction catalyzed by E. coli dihydrofolate reductase (ecDHFR). The explanation identifies two features that may be important in interpreting enzymatic kinetic isotope effects more generally. Our interpretation involves a dynamical simulation of the ecDHFR system using a combined quantum mechanical and molecular mechanical (QM/MM) approach. Special care is taken to maintain the volume consistent with a pressure of 1 atmosphere as the temperature is changed. The hydride transfer KIEs at 5°C and 45°C and 1 atm are calculated using ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT). Two features have been identified that may reduce the temperature dependence of KIEs for enzymatic reactions as compared to small-molecule reactions, namely (1) variation of the transition state position and (2) temperature dependence of the effective potential barrier for tunneling. An EA-VTST/MT calculation including these effects predicts a KIE that changes by only 6.5% over the 5–45°C interval for the hydride transfer catalyzed by ecDHFR, which is only slightly larger than the stated experimental uncertainty on the KIE. Without these two mechanisms, it would have decreased by 16%. On the one hand, our calculation is able to agree with the available experiment without introducing any new assumptions or special considerations. On the other hand our calculation provides a caution to those who invoke intrinsically temperature-independent mechanisms to discuss such reactions. We would describe the reaction as having only a small temperature dependence to the kinetic isotope effect, not as have a temperature-independent kinetic isotope effect. Since the small temperature dependence arises from competing effects, it will be extremely hard to predict it quantitatively; and the qualitative understanding of the factors affecting such KIEs is a more important take-home message than the precisely calculated value of the temperature dependence.

 

711. "Dependence of Transition State Structure on Substrate: The Intrinsic C-13 Kinetic Isotope Effect is Different for Physiological and Slow Substrate of the Ornithine Decarboxylase Reaction Because of Different Hydrogen Bonding Structures," D. Sicinska, D. G. Truhlar, and P. Paneth, Journal of the American Chemical Society 127, 5414-5422 (2005).

 

Kinetic isotope effects (KIEs) are a very sensitive tool for elucidating mechanisms of chemical and enzymatic reactions and for inferring transition state structure. Because carbon dioxide is especially amenable to convenient isotopic analysis of its carbon content, carbon-13 KIEs have been widely used in studies of spontaneous and enzymatic decarboxylation reactions. The enzymatic processes, however, are always complex reactions involving more than one step. Even when the chemical step in an enzymatic reaction is the only isotopically sensitive step, the observed kinetic isotope effect may be different from the KIE for the chemical step, which is called the intrinsic KIE. One approach to obtaining the intrinsic KIE is to use an alternative reactant that reacts much more slowly and to assume that the KIE measured for this case is the intrinsic KIE. A similar approach relies on the use of an appropriately mutated enzyme to slow the reaction. Although the slow-substrate and mutated-enzyme hypotheses have been widely invoked, sometimes justified as reasonable, and sometimes questioned, the precise limits of their validity are unknown, and examples of how they might fail (if and when they do fail) would be very instructive. Ornithine decarboxylase (ODC, EC 4.1.1.17) is a good case for testing these procedures since experimental values of carbon kinetic isotope effects on the ornithine decarboxylase reaction and several of its active-site mutants have been already reported. Ornithine decarboxylase is the first and the rate-controlling enzyme in polyamine biosynthesis; it decarboxylates L-ornithine to form the diamine putrescine. We present calculations performed using a combined quantum mechanical and molecular mechanical (QM/MM) method with the AM1 semiempirical Hamiltonian for the wild-type ornithine decarboxylase reaction with ornithine (the physiological substrate) and lysine (a "slow" substrate) and for mutant E274A with ornithine substrate. The dynamical method is variational transition state theory with quantized vibrations. We employ a single reaction coordinate equal to the carbon-carbon distance of the dissociating bond, and we find a large difference between the intrinsic kinetic isotope effect for the physiological substrate, which equals 1.04, and that for the slow substrate, which equals 1.06. This shows that, contrary to a commonly accepted assumption, kinetic isotope effects on slow substrates are not always good models of intrinsic kinetic isotope effects on physiological substrates. Furthermore, analysis of free-energy-based samples of transition state structures shows that the differences in kinetic isotope effects may be traced to different numbers of hydrogen bonds at the different transition states of the different reactions.