A generalized hybrid orbital (GHO) method has been developed at the semiempirical level in combined quantum mechanical and molecular mechanical (QM/MM) calculations. In this method, a set of hybrid orbitals is placed on the boundary atom between the QM and MM fragments, and one of the hybrid orbitals participates in the SCF calculation for the atoms in the QM region. The GHO method provides a well-defined potential energy surface for a hybrid QM/MM system and is a significant improvement over the “link-atom” approach by saturating the QM valencies with hydrogen atoms. The method has been tested on small molecules and yields reasonable structural, energetic, and electronic results in comparison with the results of the corresponding QM and MM approximations. The GHO method will greatly increase the applicability of combined QM/MM methods to systems comprising large molecules, such as proteins.

We present a theoretical framework for the calculation of rate constants of enzyme-catalyzed reactions that combines variational optimization of the dynamical bottleneck for overbarrier reactive events and multidimensional quantum mechanical tunneling dynamics for through-barrier reactive events, both in the presence of the protein environment. The theory features a two-zone, three-stage procedure called ensemble-averaged variational transition state theory with multidimensional tunneling (EA-VTST/MT) with the transmission coefficient based on the equilibrium secondary-zone (ESZ) approximation for including the effects of the protein on a catalytic reaction center, called the primary zone. The dynamics is calculated by canonical variational theory with optimized multidimensional tunneling contributions, and the formalism allows for Boltzmann averaging over an ensemble of reactant and transition state conformations. In the first stage of the calculations, we assume that the generalized transition states can be well described by a single progress coordinate expressed in primary-zone internal coordinates; in subsequent steps, the transmission coefficient is averaged over a set of primary-zone reaction paths that depend on the protein configuration, and each reaction path has its own reaction coordinate and optimized tunneling path. We also present a simpler approximation to the transmission coefficient that is called the static secondary-zone (SSZ) approximation. We illustrate both versions of this method by carrying out calculations of the reaction rate constants and kinetic isotope effects for oxidation of benzyl alcoholate to benzaldehyde by horse liver alcohol dehydrogenase. The potential energy surface is modeled by a combined generalized hybrid orbital/quantum mechanical/molecular mechanical/semiempirical valence bond (GHO-QM/MM/SEVB) method. The multidimensional tunneling calculations are microcanonically optimized by employing both the small-curvature tunneling approximation and version 4 of the large-curvature tunneling approximation. We find that the variation of the protein mean force as a function of reaction coordinate is quantitatively significant, but it does not change the qualitative conclusions for the present reaction. We obtain good agreement with experiment for both kinetic isotope effects and Swain−Schaad exponents.

One of the strongest experimental indications of hydrogen tunneling in biology has been the elevated Swain-Schaad exponent for the secondary kinetic isotope effect in the hydride-transfer step catalyzed by liver alcohol dehydrogenase. This process has been simulated using canonical variational transition-state theory for overbarrier dynamics and optimized multidimensional paths for tunneling. Semiclassical quantum effects on the dynamics are included on a 21-atom substrate-enzyme-coenzyme primary zone embedded in the potential of a substrate-enzyme-coenzyme-solvent secondary zone. The potential energy surface is calculated by treating 54 atoms by quantum mechanical electronic structure methods and 5506 protein, coenzyme, and solvent atoms by molecular mechanical force fields. We find an elevated Swain-Schaad exponent for the secondary kinetic isotope effect and generally good agreement with other experimental observables. Quantum mechanical tunneling is calculated to account for ∼60% of the reactive flux, confirming the dominance of tunneling that was inferred from the Swain-Schaad exponent. The calculations provide a detailed picture of the origin of the kinetic isotope effect and the nature of the tunneling process.