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.     

Contributions of pauli repulsions to the energetics and physical properties computed in QM/MM methods

Conventional combined quantum mechanical/molecular mechanical (QM/MM) methods lack explicit treatment of Pauli repulsions between the quantum‐mechanical and molecular‐mechanical subsystems. Instead, classical Lennard‐Jones (LJ) potentials between QM and MM nuclei are used to model electronic Pauli repulsion and long‐range London dispersion, despite the fact that the latter two are inherently of quantum nature. Use of the simple LJ potential in QM/MM methods can reproduce minimal geometries and energies of many molecular clusters reasonably well, as compared to full QM calculations. However, we show here that the LJ potential cannot correctly describe subtle details of the electron density of the QM subsystem because of the neglect of Pauli repulsions between the QM and MM subsystems. The inaccurate electron density subsequently affects the calculation of electronic and magnetic properties of the QM subsystem. To explicitly consider Pauli interactions with QM/MM methods, we propose a method to use empirical effective potentials on the MM atoms. The test case of the binding energy and magnetic properties of a water dimer shows promising results for the general application of effective potentials to mimic Pauli repulsions in QM/MM calculations. © 2013 Wiley Periodicals, Inc.     

Quantum supercharger library: Hyper‐parallel integral derivatives algorithms for ab initio QM/MM dynamics

This article describes an extension of the quantum supercharger library (QSL) to perform quantum mechanical (QM) gradient and optimization calculations as well as hybrid QM and molecular mechanical (QM/MM) molecular dynamics simulations. The integral derivatives are, after the two‐electron integrals, the most computationally expensive part of the aforementioned calculations/simulations. Algorithms are presented for accelerating the one‐ and two‐electron integral derivatives on a graphical processing unit (GPU). It is shown that a Hartree–Fock ab initio gradient calculation is up to 9.3X faster on a single GPU compared with a single central processing unit running an optimized serial version of GAMESS‐UK, which uses the efficient Schlegel method for ‐ and ‐orbitals. Benchmark QM and QM/MM molecular dynamics simulations are performed on cellobiose in vacuo and in a 39 Å water sphere (45 QM atoms and 24843 point charges, respectively) using the 6‐31G basis set. The QSL can perform 9.7 ps/day of ab initio QM dynamics and 6.4 ps/day of QM/MM dynamics on a single GPU in full double precision. © 2015 Wiley Periodicals, Inc.     

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.     

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.     

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.     

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.     

The QM‐MM interface for CHARMM‐deMon

We present a new QM/MM interface for fast and efficient simulations of organic and biological molecules. The CHARMM/deMon interface has been developed and tested to perform minimization and atomistic simulations for multi‐particle systems. The current features of this QM/MM interface include readability for molecular dynamics, tested compatibility with Free Energy Perturbation simulations (FEP) using the dual topology/single coordinate method. The current coupling scheme uses link atoms, but further extensions of the code to incorporate other available schemes are planned. We report the performance of different levels of theory for the treatment of the QM region, while the MM region was represented by a classical force‐field (CHARMM27) or a polarizable force‐field based on a simple Drude model. The current QM/MM implementation can be coupled to the dual‐thermostat method and the VV2 integrator to run molecular dynamics simulations. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Fragmentation-based QM/MM simulations: length dependence of chain dynamics and hydrogen bonding of polyethylene oxide and polyethylene in aqueous solutions

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEO(n) (n = 6-20), and polyethylene, PE(n) ( n = 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations.     

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 range of applicability of structure‐ and energy‐adjusted QM/MM link bonds

The hydrogen‐capping method is one of the most popular and widely used coupling‐schemes for quantum mechanics/molecular mechanics (QM/MM)‐molecular dynamics simulations of macromolecular systems. This is mostly due to the fact that it is fairly convenient to implement and parametrize, thus providing an excellent compromise between accuracy and computational effort. In this work, a viable and straight‐forward approach to optimize the placing of the link atom on a suitable distance ratio between the frontier atoms is discussed. To further increase the accuracy, instead of global parameters for all amino acids, different parameter sets for each type of amino acid are derived. The dependency of the link bond parameters on the chemical environment and the used QM‐method is probed to assess the range of applicability of the parametrization. Suitable sets of parameters for RI‐MP2, B3LYP, (RI)‐B3LYP‐D3, and RI‐BLYP‐D3 at triple‐zeta level for all relevant proteinogenic amino acids are presented. Furthermore, the scope and range of the perturbation, stemming from the introduction of link bonds is evaluated through application of the presented QM/MM scheme in calculations of the active site of 15S‐lipoxygenase. © 2015 Wiley Periodicals, Inc.     

Methodological aspects of QM/MM calculations: A case study on matrix metalloproteinase‐2

We address methodological issues in quantum mechanics/molecular mechanics (QM/MM) calculations on a zinc‐dependent enzyme. We focus on the first stage of peptide bond cleavage by matrix metalloproteinase‐2 (MMP‐2), that is, the nucleophilic attack of the zinc‐coordinating water molecule on the carbonyl carbon atom of the scissile fragment of the substrate. This step is accompanied by significant charge redistribution around the zinc cation, bond cleavage, and bond formation. We vary the size and initial geometry of the model system as well as the computational protocol to demonstrate the influence of these choices on the results obtained. We present QM/MM potential energy profiles for a set of snapshots randomly selected from QM/MM‐based molecular dynamics simulations and analyze the differences in the computed profiles in structural terms. Since the substrate in MMP‐2 is located on the protein surface, we investigate the influence of the thickness of the water layer around the enzyme on the QM/MM energy profile. Thin water layers (0–2 Å) give unrealistic results because of structural reorganizations in the active‐site region at the protein surface. A 12 Å water layer appears to be sufficient to capture the effect of the solvent; the corresponding QM/MM energy profile is very close to that obtained from QM/MM/SMBP calculations using the solvent macromolecular boundary potential (SMBP). We apply the optimized computational protocol to explain the origin of the different catalytic activity of the Glu116Asp mutant: the energy barrier for the first step is higher, which is rationalized on structural grounds. © 2016 Wiley Periodicals, Inc.     

Incorporation of a QM/MM buffer zone in the variational double self-consistent field method

The explicit polarization (X-Pol) potential is an electronic-structure-based polarization force field, designed for molecular dynamics simulations and modeling of biopolymers. In this approach, molecular polarization and charge transfer effects are explicitly treated by a combined quantum mechanical and molecular mechanical (QM/MM) scheme, and the wave function of the entire system is variationally optimized by a double self-consistent field (DSCF) method. In the present article, we introduce a QM buffer zone for a smooth transition from a QM region to an MM region. Instead of using the Mulliken charge approximation for all QM/MM interactions, the Coulombic interactions between the adjacent fragments are determined directly by electronic structure theory. The present method is designed to accelerate the speed of convergence of the total energy and charge density of the system.     

QM/MM investigation of the degradation mechanism of the electron-transporting layer

The mechanism of charge transfer among tris(8-hydroxyquinolinate)aluminum (Alq₃) molecules in the electron-transporting layer (ETL) under amorphous conditions was theoretically investigated using both quantum mechanical/molecular mechanical (QM/MM) calculations and molecular dynamics (MD) simulations. The rate constant of the electron transfer was estimated for the equilibrated structure taken from the QM/MM MD simulations, based on the hopping model and Marcus theory. It was found that the coordination of a (LiF)₄ cluster in ETL drastically lowers the energy of the lowest unoccupied molecular orbital in the Alq₃ molecule. The small rate constant, namely the slow charge mobility, in ETL is believed to be causally related to the low-lying delocalized unoccupied molecular orbital of Alq₃ coordinated by the (LiF)₄ cluster. The results suggest that their interaction has a considerable influence on efficiency and is attributed in part to ETL degradation in organic light-emitting diodes.     

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.     

The CHARMM–TURBOMOLE interface for efficient and accurate QM/MM molecular dynamics,free energies,and excited state properties

The quantum mechanical (QM)/molecular mechanical (MM) interface between Chemistry at HARvard Molecular Mechanics (CHARMM) and TURBOMOLE is described. CHARMM provides an extensive set of simulation algorithms, like molecular dynamics (MD) and free energy perturbation, and support for mature nonpolarizable and Drude polarizable force fields. TURBOMOLE provides fast QM calculations using density functional theory or wave function methods and excited state properties. CHARMM–TURBOMOLE is well‐suited for extended QM/MM MD simulations using first principles methods with large (triple‐ζ) basis sets. We demonstrate these capabilities with a QM/MM simulation of Mg²⁺(aq), where the MM outer sphere water molecules are represented using the SWM4‐NDP Drude polarizable force field and the ion and inner coordination sphere are represented using QM PBE, PBE0, and MP2 methods. The relative solvation free energies of Mg²⁺ and Zn²⁺ were calculated using thermodynamic integration. We also demonstrate the features for excited state properties. We calculate the time‐averaged solution absorption spectrum of indole, the emission spectrum of the indole excited state, and the electronic circular dichroism spectrum of an oxacepham. © 2014 Wiley Periodicals, Inc.     

QM/MM calculation of solvent effects on absorption spectra of guanine

Electronic spectra of guanine in the gas phase and in water were studied by quantum mechanical/molecular mechanical (QM/MM) methods. Geometries for the excited‐state calculations were extracted from ground‐state molecular dynamics (MD) simulations using the self‐consistent‐charge density functional tight binding (SCC‐DFTB) method for the QM region and the TIP3P force field for the water environment. Theoretical absorption spectra were generated from excitation energies and oscillator strengths calculated for 50 to 500 MD snapshots of guanine in the gas phase (QM) and in solution (QM/MM). The excited‐state calculations used time‐dependent density functional theory (TDDFT) and the DFT‐based multireference configuration interaction (DFT/MRCI) method of Grimme and Waletzke, in combination with two basis sets. Our investigation covered keto‐N7H and keto‐N9H guanine, with particular focus on solvent effects in the low‐energy spectrum of the keto‐N9H tautomer. When compared with the vertical excitation energies of gas‐phase guanine at the optimized DFT (B3LYP/TZVP) geometry, the maxima in the computed solution spectra are shifted by several tenths of an eV. Three effects contribute: the use of SCC‐DFTB‐based rather than B3LYP‐based geometries in the MD snapshots (red shift of ca. 0.1 eV), explicit inclusion of nuclear motion through the MD snapshots (red shift of ca. 0.1 eV), and intrinsic solvent effects (differences in the absorption maxima in the computed gas‐phase and solution spectra, typically ca. 0.1–0.3 eV). A detailed analysis of the results indicates that the intrinsic solvent effects arise both from solvent‐induced structural changes and from electrostatic solute–solvent interactions, the latter being dominant. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

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