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     

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.     

Conformational simulations of aqueous solvated alpha-conotoxin GI and its single disulfide analogues using a polarizable force field model

The solution conformation of alpha-conotoxin GI and its two single disulfide analogues are simulated using a polarizable force field in combination with the molecular fragmentation quantum chemical calculation. The polarizability is explicitly described by allowing the partial charges and fragment dipole moments to be variables, with values coming from the linear-scaling energy-based molecular fragmentation calculations at the B3LYP/6-31G(d) level. In comparison with the full quantum chemical calculations, the fragmentation approaches can yield precise ground-state energies, dipole moments, and static polarizabilities for peptides. The B3LYP/6-31G(d) charges and fragment-centered dipole moments are introduced in calculations of electrostatic terms in both AmberFF03 and OPLS force fields. Our test calculations on the gas-phase glucagon (PDB code: 1gcn) and solvated alpha-conotoxin GI (PDB code: 1not) demonstrate that the present polarization model is capable of describing the structural properties (such as the relative conformational energies, intramolecular hydrogen bonds, and disulfide bonds) with accuracy comparable to some other polarizable force fields (ABEEM/MM and OPLS-PFF) and the quantum mechanics/molecular mechanics (QM/MM) hybrid model. The employment of fragment-centered dipole moments in calculations of dipole-dipole interactions can save computational time in comparison with those polarization models using atom-centered dipole moments without much loss of accuracy. The molecular dynamics simulations using the polarizable force field demonstrate that two single disulfide GI analogues are more flexible and less structured than the native alpha-conotoxin GI, in agreement with NMR experiments. The polarization effect is important in simulations of the folding/unfolding process of solvated proteins.     

Quantum mechanical/molecular mechanical study on the mechanism of the enzymatic Baeyer-Villiger reaction

We report a combined quantum mechanical/molecular mechanical (QM/MM) study on the mechanism of the enzymatic Baeyer-Villiger reaction catalyzed by cyclohexanone monooxygenase (CHMO). In QM/MM geometry optimizations and reaction path calculations, density functional theory (B3LYP/TZVP) is used to describe the QM region consisting of the substrate (cyclohexanone), the isoalloxazine ring of C4a-peroxyflavin, the side chain of Arg-329, and the nicotinamide ring and the adjacent ribose of NADP(+), while the remainder of the enzyme is represented by the CHARMM force field. QM/MM molecular dynamics simulations and free energy calculations at the semiempirical OM3/CHARMM level employ the same QM/MM partitioning. According to the QM/MM calculations, the enzyme-reactant complex contains an anionic deprotonated C4a-peroxyflavin that is stabilized by strong hydrogen bonds with the Arg-329 residue and the NADP(+) cofactor. The CHMO-catalyzed reaction proceeds via a Criegee intermediate having pronounced anionic character. The initial addition reaction has to overcome an energy barrier of about 9 kcal/mol. The formed Criegee intermediate occupies a shallow minimum on the QM/MM potential energy surface and can undergo fragmentation to the lactone product by surmounting a second energy barrier of about 7 kcal/mol. The transition state for the latter migration step is the highest point on the QM/MM energy profile. Gas-phase reoptimizations of the QM region lead to higher barriers and confirm the crucial role of the Arg-329 residue and the NADP(+) cofactor for the catalytic efficiency of CHMO. QM/MM calculations for the CHMO-catalyzed oxidation of 4-methylcyclohexanone reproduce and rationalize the experimentally observed (S)-enantioselectivity for this substrate, which is governed by the conformational preferences of the corresponding Criegee intermediate and the subsequent transition state for the migration step.     

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.     

Probing the transition from hydrophilic to hydrophobic solvation with atomic scale resolution

Picosecond and femtosecond X-ray absorption spectroscopy is used to probe the changes of the solvent shell structure upon electron abstraction of aqueous iodide using an ultrashort laser pulse. The transient L(1,3) edge EXAFS at 50 ps time delay points to the formation of an expanded water cavity around the iodine atom, in good agreement with classical and quantum mechanical/molecular mechanics (QM/MM) molecular dynamics (MD) simulations. These also show that while the hydrogen atoms pointed toward iodide, they predominantly point toward the bulk solvent in the case of iodine, suggesting a hydrophobic behavior. This is further confirmed by quantum chemical (QC) calculations of I(-)/I(0)(H(2)O)(n=1-4) clusters. The L(1) edge sub-picosecond spectra point to the existence of a transient species that is not present at 50 ps. The QC calculations and the QM/MM MD simulations identify this transient species as an I(0)(OH(2)) complex inside the cavity. The simulations show that upon electron abstraction most of the water molecules move away from iodine, while one comes closer to form the complex that lives for 3-4 ps. This time is governed by the reorganization of the main solvation shell, basically the time it takes for the water molecules to reform an H-bond network. Only then is the interaction with the solvation shell strong enough to pull the water molecule of the complex toward the bulk solvent. Overall, much of the behavior at early times is determined by the reorientational dynamics of water molecules and the formation of a complete network of hydrogen bonded molecules in the first solvation shell.     

Characteristics of CO3(2-)-water hydrogen bonds in aqueous solution: insights from HF/MM and B3LYP/MM MD simulations

Two combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, namely HF/MM and B3LYP/MM, have been performed to investigate the local hydration structure and dynamics of carbonate (CO(3)(2-)) in dilute aqueous solution. With respect to the QM/MM scheme, the QM region, which contains the CO(3)(2-) and its surrounding water molecules, was treated at HF and B3LYP levels of accuracy, respectively, using the DZV+ basis set, while the rest of the system is described by classical MM potentials. For both the HF/MM and B3LYP/MM simulations, it is observed that the hydrogen bonds between CO(3)(2-) oxygens and their nearest-neighbor waters are relatively strong, i.e., compared to water-water hydrogen bonds in the bulk, and that the first shell of each CO(3)(2-) oxygen atom somewhat overlaps with the others, which allows migration of water molecules among the coordinating sites to exist. In addition, it is observed that first-shell waters are either "loosely" or "tightly" bound to the respective CO(3)(2-) oxygen atoms, leading to large fluctuations in the number of first-shell waters, ranging from 1 to 6 (HF/MM) and 2 to 7 (B3LYP/MM), with the prevalent value of 3. Upon comparing the HF and B3LYP methods in describing this hydrated ion, the latter is found to overestimate the hydrogen-bond strength in the CO(3)(2-)-water complexes, resulting in a slightly more compact hydration structure at each of the CO(3)(2-) oxygens.     

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.     

MM and QM/MM Modeling of Threonyl-tRNA Synthetase: Model Testing and Simulations

Aminoacyl-tRNA synthetases are centrally important enzymes in protein synthesis. We have investigated threonyl-tRNA synthetase from E. coli, complexed with reactants, using molecular mechanics and combined quantum mechanical/molecular mechanical (QM/MM) techniques. These modeling methods have the potential to provide molecular level understanding of enzyme catalytic processes. Modeling of this enzyme presents a number of challenges. The procedure of system preparation and testing is described in detail. For example, the number of metal ions at the active site, and their positions, were investigated. Molecular dynamics simulations suggest that the system is most stable when it contains only one magnesium ion, and the zinc ion is removed. Two different QM/MM methods were tested in models based on the findings of MM molecular dynamics simulations. AM1/CHARMM calculations resulted in unrealistic structures for the phosphates in this system. This is apparently due to an error of AM1. PM3/CHARMM calculations proved to be more suitable for this enzyme system. These results will provide a useful basis for future modeling investigations of the enzyme mechanism and dynamics.     

Ab initio quantum mechanical charge field (QMCF) molecular dynamics: a QM/MM – MD procedure for accurate simulations of ions and complexes

A new formalism for quantum mechanical / molecular mechanical (QM/MM) dynamics of chemical species in solution has been developed, which does not require the construction of any other potential functions except those for solvent–solvent interactions, maintains all the advantages of large simulation boxes and ensures the accuracy of ab initio quantum mechanics for all forces acting in the chemically most relevant region. Interactions between solute and more distant solvent molecules are incorporated by a dynamically adjusted force field corresponding to the actual molecular configuration of the simulated system and charges derived from the electron distribution in the solvate. The new formalism has been tested with some examples of hydrated ions, for which accurate conventional ab initio QM/MM simulations have been previously performed, and the comparison shows equivalence and in some aspects superiority of the new method. As this simulation procedure does not require any tedious construction of two-and three-body interaction potentials inherent to conventional QM/MM approaches, it opens the straightforward access to ab initio molecular dynamics simulations of any kind of solutes, such as metal complexes and other composite species in solution.     

Microhydration of Alanine in Gas Phase Studied by Quantum Chemical Method and ABEEMσπ/MM Fluctuating Charge Model

A fluctuating charge interaction potential function for alanine-water was constructed in the spirit of newly developed ABEEMσπ/MM(atom-bond electronegativity equalization method at the σπ level fused into molecular mechanics). The properties of gaseous neutral alanine-(H₂O)_n(n=1-7) clusters were systematically investigated by quantum mechanics(QM) and the constructed ABEEMσπ/MM potential, such as conformations, hydrogen bonds (H-bonds), interaction energies, charge distributions, and so on. The results of ABEEMσπ/MM model are in fair agreement with those of QM and available experimental data. For isolated alanine, compared with those of experimental structure, the average absolute deviations(AAD) of bond length and bond angle are 0.002 nm and 1.4°, respectively. For alanine-water clusters, the AAD of interaction energies and H-bond lengths are only 3.77 kJ/mol and 0.012 nm, respectively, compared to the results of MP2/aug-cc-pVDZ//MP2/6-311+G^** method. The ABEEMσπ charges fluctuate with the changing conformation of the system, and can accurately and reasonably reflect the interpolarization between water and alanine. The presented alanine-water potential function may provide a basis for further simulations on related aqueous solutions of biomolecules.     

Simulations of liquid ammonia based on the combined quantum mechanical/molecular mechanical (QM/MM) approach

Two combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics simulations, namely, HF/MM and B3LYP/MM, have been performed to investigate the local structure and dynamics of liquid ammonia. The most interesting region, a sphere containing a central reference molecule and all its nearest surrounding molecules (first coordination shell), was treated by the Hartree-Fock (HF) and hybrid density functional B3LYP methods, whereas the rest of the system was described by the classical pair potentials. On the basis of both HF and B3LYP methods, it is observed that the hydrogen bonding in this peculiar liquid is weak. The structure and dynamics of this liquid are suggested to be determined by the steric packing effects, rather than by the directional hydrogen bonding interactions. Compared to previous empirical as well as Car-Parrinello (CP) molecular dynamics studies, our QM/MM simulations provide detailed information that is in better agreement with experimental data.     

Structure of the enzyme-substrate complex for guanosine triphosphate hydrolysis by elongation factor EF-Tu: Comparison of quantum mechanics/molecular mechanics and molecular dynamics results

The equilibrium geometric configurations of the enzyme-substrate complex for guanosine triphosphate hydrolysis by elongation factor EF-Tu calculated using two theoretical approaches, a combined quantum mechanics/molecular mechanics (QM/MM) method and a molecular dynamics method, are compared. The reaction complex geometry determined by the QM/MM method is consistent with the accepted reaction mechanism, whereas, in the enzyme-substrate structure predicted by the molecular dynamics method with the CHARMM force field, the relative positions of the nucleophilic reagent (water molecules) and the base (a histidine side chain) do not correspond to the optimal reagent arrangement.     

Optimizing efficiency of perturbative Monte Carlo method

We introduce error weighting functions into the perturbative Monte Carlo method for use with a hybrid ab initio quantum mechanics/molecular mechanics (QM/MM) potential. The perturbative Monte Carlo approach introduced earlier provides a means to reduce the number of full SCF calculations in simulations using a QM/MM potential by evoking perturbation theory to calculate energy changes due to displacements of an MM molecule. The use of weighting functions, introduced here, allows an optimal number of MM molecule displacements to occur between the performance of the full self-consistent field calculations. This will allow the ab initio QM/MM approach to be applied to systems that require more accurate treatment of the QM and/or MM regions. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1632–1638, 1998     

Reliable treatment of electrostatics in combined QM/MM simulation of macromolecules

A robust approach for dealing with electrostatic interactions for spherical boundary conditions has been implemented in the QM/MM framework. The development was based on the generalized solvent boundary potential (GSBP) method proposed by Im et al. [J. Chem. Phys. 114, 2924 (2001)], and the specific implementation was applied to the self-consistent-charge density-functional tight-binding approach as the quantum mechanics (QM) level, although extension to other QM methods is straightforward. Compared to the popular stochastic boundary-condition scheme, the new protocol offers a balanced treatment between quantum mechanics/molecular mechanics (QM/MM) and MM/MM interactions; it also includes the effect of the bulk solvent and macromolecule atoms outside of the microscopic region at the Poisson-Boltzmann level. The new method was illustrated with application to the enzyme human carbonic anhydrase II and compared to stochastic boundary-condition simulations using different electrostatic treatments. The GSBP-based QM/MM simulations were most consistent with available experimental data, while conventional stochastic boundary simulations yielded various artifacts depending on different electrostatic models. The results highlight the importance of carefully treating electrostatics in QM/MM simulations of biomolecules and suggest that the commonly used truncation schemes should be avoided in QM/MM simulations, especially in simulations that involve extensive conformational samplings. The development of the GSBP-based QM/MM protocol has opened up the exciting possibility of studying chemical events in very complex biomolecular systems in a multiscale framework.