Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

QM/MM methods have been developed as a computationally feasible solution to QM simulation of chemical processes, such as enzyme-catalyzed reactions, within a more approximate MM representation of the condensed-phase environment. However, there has been no independent method for checking the quality of this representation, especially for highly nonisotropic protein environments such as those surrounding enzyme active sites. Hence, the validity of QM/MM methods is largely untested. Here we use the possibility of performing all-QM calculations at the semiempirical PM3 level with a linear-scaling method (MOZYME) to assess the performance of a QM/MM method (PM3/AMBER94 force field). Using two model pathways for the hydride-ion transfer reaction of the enzyme dihydrofolate reductase studied previously (Titmuss et al., Chem Phys Lett 2000, 320, 169-176), we have analyzed the reaction energy contributions (QM, QM/MM, and MM) from the QM/MM results and compared them with analogous-region components calculated via an energy partitioning scheme implemented into MOZYME. This analysis further divided the MOZYME components into Coulomb, resonance and exchange energy terms. For the model in which the MM coordinates are kept fixed during the reaction, we find that the MOZYME and QM/MM total energy profiles agree very well, but that there are significant differences in the energy components. Most significantly there is a large change (approximately 16 kcal/mol) in the MOZYME MM component due to polarization of the MM region surrounding the active site, and which arises mostly from MM atoms close to (<10 A) the active-site QM region, which is not modelled explicitly by our QM/MM method. However, for the model where the MM coordinates are allowed to vary during the reaction, we find large differences in the MOZYME and QM/MM total energy profiles, with a discrepancy of 52 kcal/mol between the relative reaction (product-reactant) energies. This is largely due to a difference in the MM energies of 58 kcal/mol, of which we can attribute approximately 40 kcal/mol to geometry effects in the MM region and the remainder, as before, to MM region polarization. Contrary to the fixed-geometry model, there is no correlation of the MM energy changes with distance from the QM region, nor are they contributed by only a few residues. Overall, the results suggest that merely extending the size of the QM region in the QM/MM calculation is not a universal solution to the MOZYME- and QM/MM-method differences. They also suggest that attaching physical significance to MOZYME Coulomb, resonance and exchange components is problematic. Although we conclude that it would be possible to reparameterize the QM/MM force field to reproduce MOZYME energies, a better way to account for both the effects of the protein environment and known deficiencies in semiempirical methods would be to parameterize the force field based on data from DFT or ab initio QM linear-scaling calculations. Such a force field could be used efficiently in MD simulations to calculate free energies.     

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree–Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self‐consistent field method, Møller–Plesset perturbation theory method, and time‐dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self‐consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations‐All Atom (OPLS‐AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented. © 2013 Wiley Periodicals, Inc.     

Conical intersections in solution: formulation, algorithm, and implementation with combined quantum mechanics/molecular mechanics method

The significance of conical intersections in photophysics, photochemistry, and photodissociation of polyatomic molecules in gas phase has been demonstrated by numerous experimental and theoretical studies. Optimization of conical intersections of small- and medium-size molecules in gas phase has currently become a routine optimization process, as it has been implemented in many electronic structure packages. However, optimization of conical intersections of small- and medium-size molecules in solution or macromolecules remains inefficient, even poorly defined, due to large number of degrees of freedom and costly evaluations of gradient difference and nonadiabatic coupling vectors. In this work, based on the sequential quantum mechanics and molecular mechanics (QM/MM) and QM/MM-minimum free energy path methods, we have designed two conical intersection optimization methods for small- and medium-size molecules in solution or macromolecules. The first one is sequential QM conical intersection optimization and MM minimization for potential energy surfaces; the second one is sequential QM conical intersection optimization and MM sampling for potential of mean force surfaces, i.e., free energy surfaces. In such methods, the region where electronic structures change remarkably is placed into the QM subsystem, while the rest of the system is placed into the MM subsystem; thus, dimensionalities of gradient difference and nonadiabatic coupling vectors are decreased due to the relatively small QM subsystem. Furthermore, in comparison with the concurrent optimization scheme, sequential QM conical intersection optimization and MM minimization or sampling reduce the number of evaluations of gradient difference and nonadiabatic coupling vectors because these vectors need to be calculated only when the QM subsystem moves, independent of the MM minimization or sampling. Taken together, costly evaluations of gradient difference and nonadiabatic coupling vectors in solution or macromolecules can be reduced significantly. Test optimizations of conical intersections of cyclopropanone and acetaldehyde in aqueous solution have been carried out successfully.     

Implementation of the IMOMM methodology for performing combined QM/MM molecular dynamics simulations and frequency calculations

A technique for implementing the integrated molecular orbital and molecular mechanics (IMOMM) methodology developed by Maseras and Morokuma that is used to perform combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, frequency calculations and simulations of macromolecules including explicit solvent is presented. Although the IMOMM methodology is generalized to any coordinate system, the implementation first described by Maseras and Morokuma requires that the QM and MM gradients be transformed into internal coordinates before they are added together. This coordinate transformation can be cumbersome for macromolecular systems and can become ill-defined during the course of a molecular dynamics simulation. We describe an implementation of the IMOMM method in which the QM and MM gradients are combined in the cartesian coordinate system, thereby avoiding potential problems associated with using the internal coordinate system. The implementation can be used to perform combined QM/MM molecular dynamics simulations and frequency calculations within the IMOMM framework. Finally, we have examined the applicability of thermochemical data derived from IMOMM framework. Finally, we have examined the applicability of thermochemical data derived from IMOMM frequency calculations. Received: 11 May 1998 / Accepted: 14 August 1998 / Published online: 16 November 1998     

Toward theoretical analysis of long-range proton transfer kinetics in biomolecular pumps

Motivated by the long-term goal of theoretically analyzing long-range proton transfer (PT) kinetics in biomolecular pumps, researchers made a number of technical developments in the framework of quantum mechanics-molecular mechanics (QM/MM) simulations. A set of collective reaction coordinates is proposed for characterizing the progress of long-range proton transfers; unlike previous suggestions, the new coordinates can describe PT along highly nonlinear three-dimensional pathways. Calculations using a realistic model of carbonic anhydrase demonstrated that adiabatic mapping using these collective coordinates gives reliable energetics and critical geometrical parameters as compared to minimum energy path calculations, which suggests that the new coordinates can be effectively used as reaction coordinate in potential of mean force calculations for long-range PT in complex systems. In addition, the generalized solvent boundary potential was implemented in the QM/MM framework for rectangular geometries, which is useful for studying reactions in membrane systems. The resulting protocol was found to produce water structure in the interior of aquaporin consistent with previous studies including a much larger number of explicit solvent and lipid molecules. The effect of electrostatics for PT through a membrane protein was also illustrated with a simple model channel embedded in different dielectric continuum environments. The encouraging results observed so far suggest that robust theoretical analysis of long-range PT kinetics in biomolecular pumps can soon be realized in a QM/MM framework.     

Geometry optimization based on linear response free energy with quantum mechanical/molecular mechanical method: applications to Menshutkin-type and Claisen rearrangement reactions in aqueous solution

The authors present a method based on a linear response theory that allows one to optimize the geometries of quantum mechanical/molecular mechanical (QM/MM) systems on the free energy surfaces. Two different forms of linear response free energy functionals are introduced, and electronic wave functions of the QM region, as well as the responses of electrostatic and Lennard-Jones potentials between QM and MM regions, are self-consistently determined. The covariant matrix relating the QM charge distribution to the MM response is evaluated by molecular dynamics (MD) simulation of the MM system. The free energy gradients with respect to the QM atomic coordinates are also calculated using the MD trajectory results. They apply the present method to calculate the free energy profiles of Menshutkin-type reaction of NH3 with CH3Cl and Claisen rearrangement of allyl vinyl ether in aqueous solution. For the Menshutkin reaction, the free energy profile calculated with the modified linear response free energy functional is in good agreement with that by the free energy perturbation calculations. They examine the nonequilibrium solvation effect on the transmission coefficient and the kinetic isotope effect for the Claisen rearrangement.     

Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach

The extent to which accuracy of electric charges plays a role in protein-ligand docking is investigated through development of a docking algorithm, which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations. In this algorithm, fixed charges of ligands obtained from force field parameterization are replaced by QM/MM calculations in the protein environment, treating only the ligands as the quantum region. The algorithm is tested on a set of 40 cocrystallized structures taken from the Protein Data Bank (PDB) and provides strong evidence that use of nonfixed charges is important. An algorithm, dubbed "Survival of the Fittest" (SOF) algorithm, is implemented to incorporate QM/MM charge calculations without any prior knowledge of native structures of the complexes. Using an iterative protocol, this algorithm is able in many cases to converge to a nativelike structure in systems where redocking of the ligand using a standard fixed charge force field exhibits nontrivial errors. The results demonstrate that polarization effects can play a significant role in determining the structures of protein-ligand complexes, and provide a promising start towards the development of more accurate docking methods for lead optimization applications.     

Convergence in the QM‐only and QM/MM modeling of enzymatic reactions: A case study for acetylene hydratase

We report systematic quantum mechanics‐only (QM‐only) and QM/molecular mechanics (MM) calculations on an enzyme‐catalyzed reaction to assess the convergence behavior of QM‐only and QM/MM energies with respect to the size of the chosen QM region. The QM and MM parts are described by density functional theory (typically B3LYP/def2‐SVP) and the CHARMM force field, respectively. Extending our previous work on acetylene hydratase with QM regions up to 157 atoms (Liao and Thiel, J. Chem. Theory Comput. 2012, 8, 3793), we performed QM/MM geometry optimizations with a QM region M4 composed of 408 atoms, as well as further QM/MM single‐point calculations with even larger QM regions up to 657 atoms. A charge deletion analysis was conducted for the previously used QM/MM model ( M3a , with a QM region of 157 atoms) to identify all MM residues with strong electrostatic contributions to the reaction energetics (typically more than 2 kcal/mol), which were then included in M4 . QM/MM calculations with this large QM region M4 lead to the same overall mechanism as the previous QM/MM calculations with M3a , but there are some variations in the relative energies of the stationary points, with a mean absolute deviation (MAD) of 2.7 kcal/mol. The energies of the two relevant transition states are close to each other at all levels applied (typically within 2 kcal/mol), with the first (second) one being rate‐limiting in the QM/MM calculations with M3a ( M4 ). QM‐only gas‐phase calculations give a very similar energy profile for QM region M4 (MAD of 1.7 kcal/mol), contrary to the situation for M3a where we had previously found significant discrepancies between the QM‐only and QM/MM results (MAD of 7.9 kcal/mol). Extension of the QM region beyond M4 up to M7 (657 atoms) leads to only rather small variations in the relative energies from single‐point QM‐only and QM/MM calculations (MAD typically about 1–2 kcal/mol). In the case of acetylene hydratase, a model with 408 QM atoms thus seems sufficient to achieve convergence in the computed relative energies to within 1–2 kcal/mol.Copyright © 2013 Wiley Periodicals, Inc.     

QM/MM modeling of enantioselective pybox-ruthenium- and box-copper-catalyzed cyclopropanation reactions: scope, performance, and applications to ligand design

An extensive comparison of full-QM (B3LYP) and QM/MM (B3LYP:UFF) levels of theory has been made for two enantioselective catalytic systems, namely, Pybox-Ru and Box-Cu complexes, in the cyclopropanation of alkenes (ethylene and styrene) with methyl diazoacetate. The geometries of the key reaction intermediates and transition structures calculated at the QM/MM level are generally in satisfactory agreement with full-QM calculated geometries. More importantly, the relative energies calculated at the QM/MM level are in good agreement with those calculated at the full-QM level in all cases. Furthermore, the QM/MM energies are often in better agreement with the stereoselectivity experimentally observed, and this suggests that QM/MM calculations can be superior to full-QM calculations when subtle differences in inter- and intramolecular interactions are important in determining the selectivity, as is the case in enantioselective catalysis. The predictive value of the model presented is validated by the explanation of the unusual enantioselectivity behavior exhibited by a new bis-oxazoline ligand, the stereogenic centers of which are quaternary carbon atoms.     

The adaptive buffered force QM/MM method in the CP2K and AMBER software packages

The implementation and validation of the adaptive buffered force (AdBF) quantum‐mechanics/molecular‐mechanics (QM/MM) method in two popular packages, CP2K and AMBER are presented. The implementations build on the existing QM/MM functionality in each code, extending it to allow for redefinition of the QM and MM regions during the simulation and reducing QM‐MM interface errors by discarding forces near the boundary according to the buffered force‐mixing approach. New adaptive thermostats, needed by force‐mixing methods, are also implemented. Different variants of the method are benchmarked by simulating the structure of bulk water, water autoprotolysis in the presence of zinc and dimethyl‐phosphate hydrolysis using various semiempirical Hamiltonians and density functional theory as the QM model. It is shown that with suitable parameters, based on force convergence tests, the AdBF QM/MM scheme can provide an accurate approximation of the structure in the dynamical QM region matching the corresponding fully QM simulations, as well as reproducing the correct energetics in all cases. Adaptive unbuffered force‐mixing and adaptive conventional QM/MM methods also provide reasonable results for some systems, but are more likely to suffer from instabilities and inaccuracies. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

QM/MM vibrational mode tracking

Vibrational spectroscopy is a powerful tool to investigate the structure and dynamics of biomolecules. When small subsystems of large molecules such as active centers of enzymes are studied, quantum chemical calculations based on quantum mechanics/molecular mechanics (QM/MM) coupling schemes are a valuable means to interpret the spectra. The goal of this work is a methodological pilot study on how to selectively and thus efficiently extract certain vibrational information for extended molecular systems described by QM/MM methods. This is achieved by an extension of the mode tracking algorithm and a comparison with the partial Hessian diagonalization approach. After validating the methodology for the CO stretching vibration of 2-butanone and a delocalized CO stretch in acetylacetone, the stretching and bending modes of the CO ligand in CO myoglobin are tracked. Such systems represent an ideal application for mode tracking, because only a few strongly localized vibrations are sought for, while the large remainder of the molecule is of interest only as far as it affects these local vibrations. This influence is treated exactly by mode tracking.     

Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide

Methodology is discussed for mixed ab initio quantum mechanics/molecular mechanics modeling of systems where the quantum mechanics (QM) and molecular mechanics (MM) regions are within the same molecule. The ab initio QM calculations are at the restricted Hartree–Fock level using the pseudospectral method of the Jaguar program while the MM part is treated with the OPLS force fields implemented in the IMPACT program. The interface between the QM and MM regions, in particular, is elaborated upon, as it is dealt with by “breaking” bonds at the boundaries and using Boys-localized orbitals found from model molecules in place of the bonds. These orbitals are kept frozen during QM calculations. Results from tests of the method to find relative conformational energies and geometries of alanine dipeptides and alanine tetrapeptides are presented along with comparisons to pure QM and pure MM calculations. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1468–1494, 1999     

QMX: A versatile environment for hybrid calculations applied to the grafting of Al2Cl3Me3 on a silica surface

We present a new software to easily perform QM:MM and QM:QM' calculations called QMX. It follows the subtraction scheme and it is implemented in the Atomic Simulation Environment (ASE). Special attention is paid to couple molecular calculations with periodic boundaries approaches. QMX inherits the flexibility and versatility of the ASE package: any combination of methods namely force field, semiempirical, first principle, and ab initio, can be used as hybrid potential energy surface (PES). Its ease of use is demonstrated by considering the adsorption of Al₂Cl₃Me₃ on silica surface and by combining different levels of theory (from standard DFT to MP2 calculations) for the so‐called High Level cluster with standard PW91 density functional theory calculations for the Low Level environment. It is shown that the High Level cluster must contain the silanol group close to the aluminum atoms. The bridging adsorption is favored by 58 kJ mol^?1 at the MP2:PW91 level with respect to the terminal position. Using large clusters at the MP2:PW91 level, it is shown that PW91 calculations are sufficient for structure optimization but that embedded methods are required for accurate energy profiles. © 2013 Wiley Periodicals, Inc.