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