Optimal scaling factors for CM1 and CM3 atomic charges in RM1-based aqueous simulations

Scaling factors for atomic charges derived from the RM1 semiempirical quantum mechanical wavefunction in conjunction with CM1 and CM3 charge models have been optimized by minimizing errors in absolute free energies of hydration, ΔG_hyd, for a set of 40 molecules. Monte Carlo statistical mechanics simulations and free energy perturbation theory were used to annihilate the solutes in gas and in a box of TIP4P water molecules. Lennard–Jones parameters from the optimized potentials for liquid simulations‐all atom (OPLS–AA) force field were utilized for the organic compounds. Optimal charge scaling factors have been determined as 1.11 and 1.14 for the CM1R and CM3R methods, respectively, and the corresponding unsigned average errors in ΔG_hyd relative to experiment were 2.05 and 1.89 kcal/mol. Computed errors in aniline and two derivatives were particularly large for RM1 and their removal from the data set lowered the overall errors to 1.61 and 1.75 kcal/mol for CM1R and CM3R. Comparisons are made to the AM1 method which yielded total errors in ΔG_hyd of 1.50 and 1.64 kcal/mol for CM1A*1.14 and CM3A*1.15, respectively. This work is motivated by the need for a highly efficient yet accurate quantum mechanical (QM) method to study condensed‐phase and enzymatic chemical reactions via mixed QM and molecular mechanical (QM/MM) simulations. As an initial test, the Menshutkin reaction between NH₃ and CH₃Cl in water was computed using a RM1/TIP4P‐Ew/CM3R procedure and the resultant ΔG^?, ΔG_rxn, and geometries were in reasonable accord with other computational methods; however, some potentially serious shortcomings in RM1 are discussed. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Mechanism of proteolysis in matrix metalloproteinase‐2 revealed by QM/MM modeling

The mechanism of enzymatic peptide hydrolysis in matrix metalloproteinase‐2 (MMP‐2) was studied at atomic resolution through quantum mechanics/molecular mechanics (QM/MM) simulations. An all‐atom three‐dimensional molecular model was constructed on the basis of a crystal structure from the Protein Data Bank (ID: 1QIB), and the oligopeptide Ace‐Gln‐Gly～Ile‐Ala‐Gly‐Nme was considered as the substrate. Two QM/MM software packages and several computational protocols were employed to calculate QM/MM energy profiles for a four‐step mechanism involving an initial nucleophilic attack followed by hydrogen bond rearrangement, proton transfer, and C? N bond cleavage. These QM/MM calculations consistently yield rather low overall barriers for the chemical steps, in the range of 5–10 kcal/mol, for diverse QM treatments (PBE0, B3LYP, and BB1K density functionals as well as local coupled cluster treatments) and two MM force fields (CHARMM and AMBER). It, thus, seems likely that product release is the rate‐limiting step in MMP‐2 catalysis. This is supported by an exploration of various release channels through QM/MM reaction path calculations and steered molecular dynamics simulations. © 2015 Wiley Periodicals, Inc.     

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.     

Calculation of wave‐functions with frozen orbitals in mixed quantum mechanics/molecular mechanics methods. II. Application of the local basis equation

The application of the local basis equation (Ferenczy and Adams, J. Chem. Phys. 2009 , 130, 134108) in mixed quantum mechanics/molecular mechanics (QM/MM) and quantum mechanics/quantum mechanics (QM/QM) methods is investigated. This equation is suitable to derive local basis nonorthogonal orbitals that minimize the energy of the system and it exhibits good convergence properties in a self‐consistent field solution. These features make the equation appropriate to be used in mixed QM/MM and QM/QM methods to optimize orbitals in the field of frozen localized orbitals connecting the subsystems. Calculations performed for several properties in divers systems show that the method is robust with various choices of the frozen orbitals and frontier atom properties. With appropriate basis set assignment, it gives results equivalent with those of a related approach [G. G. Ferenczy previous paper in this issue] using the Huzinaga equation. Thus, the local basis equation can be used in mixed QM/MM methods with small size quantum subsystems to calculate properties in good agreement with reference Hartree–Fock–Roothaan results. It is shown that bond charges are not necessary when the local basis equation is applied, although they are required for the self‐consistent field solution of the Huzinaga equation based method. Conversely, the deformation of the wave‐function near to the boundary is observed without bond charges and this has a significant effect on deprotonation energies but a less pronounced effect when the total charge of the system is conserved. The local basis equation can also be used to define a two layer quantum system with nonorthogonal localized orbitals surrounding the central delocalized quantum subsystem. © 2013 Wiley Periodicals, Inc.     

CPMD/GULP QM/MM interface for modeling periodic solids: Implementation and its application in the study of Y‐zeolite supported Rhn clusters

We report here the development of hybrid quantum mechanics/molecular mechanics (QM/MM) interface between the plane‐wave density functional theory based CPMD code and the empirical force‐field based GULP code for modeling periodic solids and surfaces. The hybrid QM/MM interface is based on the electrostatic coupling between QM and MM regions. The interface is designed for carrying out full relaxation of all the QM and MM atoms during geometry optimizations and molecular dynamics simulations, including the boundary atoms. Both Born–Oppenheimer and Car–Parrinello molecular dynamics schemes are enabled for the QM part during the QM/MM calculations. This interface has the advantage of parallelization of both the programs such that the QM and MM force evaluations can be carried out in parallel to model large systems. The interface program is first validated for total energy conservation and parallel scaling performance is benchmarked. Oxygen vacancy in α‐cristobalite is then studied in detail and the results are compared with a fully QM calculation and experimental data. Subsequently, we use our implementation to investigate the structure of rhodium cluster (Rh_n; n = 2 to 6) formed from Rh(C₂H₄)₂ complex adsorbed within a cavity of Y‐zeolite in a reducible atmosphere of H₂ gas. © 2016 Wiley Periodicals, Inc.     

Quantum mechanically derived AMBER-compatible heme parameters for various states of the cytochrome P450 catalytic cycle

Molecular mechanics (MM) methods are computationally affordable tools for screening chemical libraries of novel compounds for sites of P450 metabolism. One challenge for MM methods has been the absence of a consistent and transferable set of parameters for the heme within the P450 active site. Experimental data indicate that mammalian P450 enzymes vary greatly in the size, architecture, and plasticity of their active sites. Thus, obtaining X-ray-based geometries for the development of accurate MM parameters for the major classes of hepatic P450 remains a daunting task. Our previous work with preliminary gas-phase quantum mechanics (QM)-derived atomic partial charges greatly improved the accuracy of docking studies of raloxifene to CYP3A4. We have therefore developed and tested a consistent set of transferable MM parameters based on gas-phase QM calculations of two model systems of the heme-a truncated (T-HM) and a full (F-HM) for four states of the P450 catalytic cycle. Our results indicate that the use of the atomic partial charges from the F-HM further improves the accuracy of docked predictions for raloxifene to CYP3A4. Different patterns for substrate docking are also observed depending on the choice of heme model and state. Newly parameterized heme models are tested in implicit and explicitly solvated MD simulations in the absence and presence of enzyme structures, for CYP3A4, and appear to be stable on the nanosecond simulation timescale. The new force field for the various heme states may aid the community for simulations of P450 enzymes and other heme-containing enzymes.     

QM:QM electronic embedding using Mulliken atomic charges: energies and analytic gradients in an ONIOM framework

An accurate first-principles treatment of chemical reactions for large systems remains a significant challenge facing electronic structure theory. Hybrid models, such as quantum mechanics:molecular mechanics (QM:MM) and quantum mechanics:quantum mechanics (QM:QM) schemes, provide a promising avenue for such studies. For many chemistries, including important reactions in materials science, molecular mechanics or semiempirical methods may not be appropriate, or parameters may not be available (e.g., surface chemistry of compound semiconductors such as indium phosphide or catalytic chemistry of transition metal oxides). In such cases, QM:QM schemes are of particular interest. In this work, a QM:QM electronic embedding model within the ONIOM (our own N-layer integrated molecular orbital molecular mechanics) extrapolation framework is presented. To define the embedding potential, we choose the real-system low-level Mulliken atomic charges. This results in a set of well-defined and unique embedding charges. However, the parametric dependence of the charges on molecular geometry complicates the energy gradient that is necessary for the efficient exploration of potential energy surfaces. We derive an efficient form for the forces where a single set of self-consistent field response equations is solved. Initial tests of the method and key algorithmic issues are discussed.     

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.     

Solvents effects on the mechanism of cellulose hydrolysis: A QM/MM study

This article reports a combined quantum mechanics/molecular mechanics (QM/MM) investigation on the acid hydrolysis of cellulose in water using two different models, cellobiose and a 40‐unit cellulose chain. The explicitly treated solvent molecules strongly influence the conformations, intramolecular hydrogen bonds, and exoanomeric effects in these models. As these features are largely responsible for the barrier to cellulose hydrolysis, the present QM/MM results for the pathways and reaction intermediates in water are expected to be more realistic than those from a former density functional theory (DFT) study with implicit solvent (CPCM). However, in a qualitative sense, there is reasonable agreement between the DFT/CPCM and QM/MM predictions for the reaction mechanism. Differences arise mainly from specific solute–solvent hydrogen bonds that are only captured by QM/MM and not by DFT/CPCM. © 2015 Wiley Periodicals, Inc.     

Converging ligand‐binding free energies obtained with free‐energy perturbations at the quantum mechanical level

In this article, the convergence of quantum mechanical (QM) free‐energy simulations based on molecular dynamics simulations at the molecular mechanics (MM) level has been investigated. We have estimated relative free energies for the binding of nine cyclic carboxylate ligands to the octa‐acid deep‐cavity host, including the host, the ligand, and all water molecules within 4.5 Å of the ligand in the QM calculations (158–224 atoms). We use single‐step exponential averaging (ssEA) and the non‐Boltzmann Bennett acceptance ratio (NBB) methods to estimate QM/MM free energy with the semi‐empirical PM6‐DH2X method, both based on interaction energies. We show that ssEA with cumulant expansion gives a better convergence and uses half as many QM calculations as NBB, although the two methods give consistent results. With 720,000 QM calculations per transformation, QM/MM free‐energy estimates with a precision of 1 kJ/mol can be obtained for all eight relative energies with ssEA, showing that this approach can be used to calculate converged QM/MM binding free energies for realistic systems and large QM partitions. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

How Photoisomerization Drives Peptide Folding and Unfolding: Insights from QM/MM and MM Dynamics Simulations

Photoswitchable azobenzene cross‐linkers can control the folding and unfolding of peptides by photoisomerization and can thus regulate peptide affinities and enzyme activities. Using quantum mechanics/molecular mechanics (QM/MM) methods and classical MM force fields, we report the first molecular dynamics simulations of the photoinduced folding and unfolding processes in the azobenzene cross‐linked FK‐11 peptide. We find that the interactions between the peptide and the azobenzene cross‐linker are crucial for controlling the evolution of the secondary structure of the peptide and responsible for accelerating the folding and unfolding events. They also modify the photoisomerization mechanism of the azobenzene cross‐linker compared with the situation in vacuo or in solution.     

A numerically stable restrained electrostatic potential charge fitting method

Inspired by the idea of charge decomposition in calculation of the dipole preserving and polarization consistent charges (Zhang et al., J. Comput. Chem. 2011, 32, 2127), we have proposed a numerically stable restrained electrostatic potential (ESP)‐based charge fitting method for protein. The atomic charge is composed of two parts. The dominant part is fixed to a predefined value (e.g., AMBER charge), and the residual part is to be determined by restrained fitting to residual ESP on grid points around the molecule. Nonuniform weighting factors as a function of the dominant charge are assigned to the atoms. Because the residual part is several folds to several orders smaller than the dominant part, the impact of ill‐conditioning is alleviated. This charge fitting method can be used in quantum mechanical/molecular mechanical (QM/MM) simulations and similar studies, where QM calculated electronic properties are frequently mapped to partial atomic charges. © 2012 Wiley Periodicals, Inc.     

Flexible-boundary QM/MM calculations: II. Partial charge transfer across the QM/MM boundary that passes through a covalent bond

Recently, based on the principle of electronic chemical potential equalization and the principle of charge conservation, we proposed a flexible-boundary scheme that allows both partial charge transfer and self-consistent polarization between the quantum mechanical (QM) and molecular mechanical (MM) subsystems in QM/MM calculations; the scheme was applied to study the atomic charges in selected ion–solvent complexes. In the present contribution, we further extend the flexible-boundary treatment to handle the QM/MM boundary passing through covalent bonds. We find that the flexible-boundary redistributed charge and dipole schemes yield reasonable agreement with full-QM calculations for a number of molecular ions and amino acids with charged side chains. Using the full-QM results as reference, the mean unsigned deviations are computed to be 0.06 e for atomic partial charges of the QM atoms, 0.11 e for the amounts of charge transfer between the QM and MM subsystems, and 0.016 Å for the lengths of the covalent bonds that directly connect the QM and MM subsystems. The results indicate the importance of accounting for partial charge transfer across the QM/MM boundary when the QM subsystems are charged.     

Influence of electronic effect on catalytic activity of salicylaldiminato nickel(II) complexes

A means of correlation between the activity and the net charge, the metal atom net charge correlation (MANCC), which was successful in the activity prediction of the early‐transition metal catalysts, has been used to study the catalytic activities of salicylaldiminato Ni(II) complexes, the late‐transition metal catalysts, in olefin polymerization or oligomerization. A comparison with the available quantum mechanics/molecular mechanics (QM/MM) calculation data suggests that even without a detailed mechanism, MANCC results mostly agree with QM/MM calculation data regarding insertion barrier data and enthalpy change. Eight experimental complexes were further built up by modeling; their catalytic activities predominantly increased in line with the net charges on the metal atoms. The same results were obtained for the other four complexes synthesized in the present work. © 2004 Wiley Periodicals, Inc. J Polym Sci Part A: Polym Chem 42: 4765–4774, 2004     

Smooth solvation method for d-orbital semiempirical calculations of biological reactions. 2. Application to transphosphorylation thio effects in solution

Density-functional and semiempirical quantum methods and continuum dielectric and explicit solvation models are applied to study the role of solvation on the stabilization of native and thio-substituted transphosphorylation reactions. Extensive comparison is made between results obtained from the different methods. For the semiempirical methods, explicit solvation was treated using a hybrid quantum mechanical/molecular mechanical (QM/MM) approach and the implicit solvation was treated using a recently developed smooth solvation model implemented into a d-orbital semiempirical framework (MNDO/d-SCOSMO) within CHARMM. The different quantum and solvation methods were applied to the transesterification of a 3'-ribose,5'-methyl phosphodiester that serves as a nonenzymatic model for the self-cleavage reaction catalyzed by the hammerhead and hairpin ribozymes. Thio effects were studied for a double sulfur substitution at the nonbridging phosphoryl oxygen positions. The reaction profiles of both the native and double sulfur-substituted reactions from the MNDO/d-SCOSMO calculations were similar to the QM/MM results and consistent with the experimentally observed trends. These results underscore the need for a d-orbital semiempirical representation for phosphorus and sulfur for the study of experimentally observed thio effects in enzymatic and nonenzymatic phosphoryl transfer reactions. One of the major advantages of the present approach is that it can be applied to model chemical reactions at a significantly lower computational cost than either the density-functional calculations with implicit solvation or the semiempirical QM/MM simulations with explicit solvent.     

An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.