Comment on “a minimal implementation of the AMBER–GAUSSIAN interface for Ab Initio QM/MM‐MD simulation”

We comment upon the recent critique of use of the Program for User Package Interfacing and Linking (PUPIL) system for linking AMBER and GAUSSIAN in a multiscale quantum mechanical/molecular mechanics (QM/MM) simulation (Okamoto et al., J. Comput. Chem. 2011 , 32, 932). Specifically, their method for computing forces on the MM particles from the QM region via the GAUSSIAN‐03 electrical field was already implemented in PUPIL version 1.3, publicly available beginning December 2009. Some other doubtful characterizations of PUPIL are discussed briefly in the context of current awareness of open‐source codes more generally. © 2012 Wiley Periodicals, Inc.     

Incorporation of deMon2k as a new parallel quantum mechanical code for the PUPIL system

The PUPIL system is a combination of software and protocols for the systematic linkage and interoperation of molecular dynamics and quantum mechanics codes to perform QM/MD (sometimes called QM/MM) calculations. The Gaussian03 and Amber packages were added to the PUPIL suite recently. However, efficient parallel QM codes are critical because calculation of the QM forces is the overwhelming majority of the computational load. Here we report details of incorporation of the deMon2k density functional suite as a new parallel QM code. An additional motivation is to add a highly optimized, purely DFT code. We illustrate with a demonstration study of the influence of perchlorate as a dopant ion of the poly(3,4‐ethylenedioxythiophene) conducting polymer in explicit acetonitrile solvent using Amber and deMon2k. We discuss unanticipated requirements for use of a scheme for semi‐empirical correction of Kohn‐Sham eigenvalues to give physically meaningful one‐electron gap energies. We provide comparison of both geometric parameters and electronic properties for nondoped and doped systems. We also present results comparing deMon2k and Gaussian03 calculation of forces for a short sequence of steps. We discuss briefly some difficult problems of quantum zone SCF convergence for the anionically doped system. The difficulties seem to be caused by well‐know deficiencies in simple approximate exchange‐correlation functionals. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

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.     

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.     

Hybrid molecular dynamics-quantum mechanics simulations of solute spectral properties in the condensed phase: evaluation of simulation parameters

We have explored the impact of a number of basic simulation parameters on the results of a recently developed hybrid molecular dynamics-quantum mechanics (MD-QM) method (Mercer et al., J Phys Chem B 1999, 103, 7720). The method utilizes MD simulations to explore the ground-state configuration space of the system and QM evaluation of those structures to yield the time-dependent electronic transition energy, which is transformed into the optical line-broadening function using the second-order cumulant expansion. Both linear and nonlinear optical spectra can then be generated for comparison to experiment. The dependence of the resulting spectra on the length of the MD trajectory, the QM sampling rate, and the QM model chemistry have all been examined. In particular, for the system of oxazine-4 in methanol studied here, at least 20 ps of MD trajectory are needed for qualitative convergence of linear spectral properties, and >100 ps is needed for quantitative convergence. Surprisingly, little difference is found between the 3-21G and 6-31G(d) basis sets, and the CIS and TD-B3LYP methods yield remarkably similar spectra. The semiempirical INDO/s method yields the most accurate results, reproducing the experimental Stokes shift to within 5% and the FWHM to within 20%. Nonlinear 3-pulse photon echo peak shift (3PEPS) decays have also been simulated. Decays are generally poorly reproduced, though the initial peak shift which depends on the overall coupling of motions to the solute transition energy is within 15% of experiment for all model chemistries other than those using the STO-3G basis.     

A minimal implementation of the AMBER-GAUSSIAN interface for ab initio QM/MM-MD simulation

For applying to a number of theoretical methodologies based on an ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics method connecting AMBER9 with GAUSSIAN03, we have developed an AMBER-GAUSSIAN interface (AG-IF), which can be one of the simplest architectures. In the AG-IF, only a few subroutines addition is necessary to retrieve the QM/MM energy and forces, obtained by GAUSSIAN, for solving a set of Newtonian equations of motion in AMBER. It is, therefore, easy to be modified for individual applications since AG-IF utilizes most of those functions originally equipped not only in AMBER but also in GAUSSIAN. In the present minimal implementation, only AMBER is modified, whereas GAUSSIAN is left unchanged. Moreover, a different method of calculating electrostatic forces of MM atoms interacting with QM region is proposed. Using the AG-IF, we also demonstrate three examples of application: (i) the QM versus MM comparison in the radial distribution function, (ii) the free energy gradient method, and (iii) the charge from interaction energy and forces.     

molUP: A VMD plugin to handle QM and ONIOM calculations using the gaussian software

The notable advances obtained by computational (bio)chemistry provided its widespread use in many areas of science, in particular, in the study of reaction mechanisms. These studies involve a huge number of complex calculations, which are often carried out using the Gaussian suite of programs. The preparation of input files and the analysis of the output files are not easy tasks and often involve laborious and complex steps. Taking this into account, we developed molUP: a VMD plugin that offers a complete set of tools that enhance the preparation of QM and ONIOM (QM/MM, QM/QM, and QM/QM/MM) calculations. The starting structures for these calculations can be imported from different chemical formats. A set of tools is available to help the user to examine or modify any geometry parameter. This includes the definition of layers in ONIOM calculations, choosing fixed atoms during geometry optimizations, the recalculation or adjustment of the atomic charges, performing SCANs or IRC calculations, etc. molUP also extracts the geometries from the output files as well as the energies of each of them. All of these tasks are performed in an interactive GUI that is extremely helpful for the user. MolUP was developed to be easy to handle by inexperienced users, but simultaneously to be a fast and flexible graphical interface to allow the advanced users to take full advantage of this plugin. The program is available, free of charges, for macOS, Linux, and Windows at the PortoBioComp page https://www.fc.up.pt/PortoBioComp/database/doku.php?id=molup . © 2018 Wiley Periodicals, Inc.     

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.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Towards the Accurate Thermodynamic Characterization of Enzyme Reaction Mechanisms

We employed QM/MM molecular dynamics (MD) simulations to characterize the rate-limiting step of the glycosylation reaction of pancreatic α-amylase with combined DFT/molecular dynamics methods (PBE/def2-SVP : AMBER). Upon careful choice of four starting active site conformations based on thorough reactivity criteria, Gibbs energy profiles were calculated with umbrella sampling simulations within a statistical convergence of 1–2 kcal ⋅ mol⁻¹. Nevertheless, Gibbs activation barriers and reaction energies still varied from 11.0 to 16.8 kcal ⋅ mol⁻¹ and −6.3 to +3.8 kcal ⋅ mol⁻¹ depending on the starting conformations, showing that despite significant state-of-the-art QM/MM MD sampling (0.5 ns/profile) the result still depends on the starting structure. The results supported the one step dissociative mechanism of Asp197 glycosylation preceded by an acid-base reaction by the Glu233, which are qualitatively similar to those from multi-PES QM/MM studies, and thus support the use of the latter to determine enzyme reaction mechanisms.     

Flooding in GROMACS: accelerated barrier crossings in molecular dynamics

The major bottleneck of today's atomistic molecular dynamics (MD) simulations is that because of the enormous computational effort involved, only processes at nanoseconds to microseconds time scales or faster can be studied directly. Unfortunately, apart from a few exceptions, relevant processes, such as chemical reactions or many large scale conformational transitions in proteins, occur at slower time scales and therefore are currently far out of reach for conventional MD. The flooding technique addresses this problem by inclusion of a flooding potential into the force field. This flooding potential locally destabilizes the educt state and thereby significantly accelerates the escape from the initial energy well without affecting the reaction pathway. Here, we summarize the theory and method for the computational chemistry community and detail the implementation within the official version 3.3 of the freely available MD program package GROMACS. Two examples shall demonstrate the application of flooding to accelerate conformational transitions and chemical reactions. The second example was carried out within a QM/MM framework.     

Computer-aided drug design platform using PyMOL

The understanding and optimization of protein-ligand interactions are instrumental to medicinal chemists investigating potential drug candidates. Over the past couple of decades, many powerful standalone tools for computer-aided drug discovery have been developed in academia providing insight into protein-ligand interactions. As programs are developed by various research groups, a consistent user-friendly graphical working environment combining computational techniques such as docking, scoring, molecular dynamics simulations, and free energy calculations is needed. Utilizing PyMOL we have developed such a graphical user interface incorporating individual academic packages designed for protein preparation (AMBER package and Reduce), molecular mechanics applications (AMBER package), and docking and scoring (AutoDock Vina and SLIDE). In addition to amassing several computational tools under one interface, the computational platform also provides a user-friendly combination of different programs. For example, utilizing a molecular dynamics (MD) simulation performed with AMBER as input for ensemble docking with AutoDock Vina. The overarching goal of this work was to provide a computational platform that facilitates medicinal chemists, many who are not experts in computational methodologies, to utilize several common computational techniques germane to drug discovery. Furthermore, our software is open source and is aimed to initiate collaborative efforts among computational researchers to combine other open source computational methods under a single, easily understandable graphical user interface.     

Reaction mechanism of Ru(II) piano‐stool complexes: Umbrella sampling QM/MM MD study

Biologically relevant interactions of piano‐stool ruthenium(II) complexes with ds‐DNA are studied in this article by hybrid quantum mechanics—molecular mechanics (QM/MM) computational technique. The whole reaction mechanism is divided into three phases: (i) hydration of the [Ru^II(η⁶‐benzene)(en)Cl]⁺ complex, (ii) monoadduct formation between the resulting aqua‐Ru(II) complex and N7 position of one of the guanines in the ds‐DNA oligomer, and (iii) formation of the intrastrand Ru(II) bridge (cross‐link) between two adjacent guanines. Free energy profiles of all the reactions are explored by QM/MM MD umbrella sampling approach where the Ru(II) complex and two guanines represent a quantum core, which is described by density functional theory methods. The combined QM/MM scheme is realized by our own software, which was developed to couple several quantum chemical programs (in this study Gaussian 09) and Amber 11 package. Calculated free energy barriers of the both ruthenium hydration and Ru(II)‐N7(G) DNA binding process are in good agreement with experimentally measured rate constants. Then, this method was used to study the possibility of cross‐link formation. One feasible pathway leading to Ru(II) guanine‐guanine cross‐link with synchronous releasing of the benzene ligand is predicted. The cross‐linking is an exergonic process with the energy barrier lower than for the monoadduct reaction of Ru(II) complex with ds‐DNA. © 2014 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.     

Computational characterization of reaction intermediates in the photocycle of the sensory domain of the AppA blue light photoreceptor

The AppA protein with the BLUF (blue light using flavin adenine dinucleotide) domain is a blue light photoreceptor that cycle between dark-adapted and light-induced functional states. We characterized possible reaction intermediates in the photocycle of AppA BLUF. Molecular dynamics (MD), quantum chemical and quantum mechanical-molecular mechanical (QM/MM) calculations were carried out to describe several stable structures of a molecular system modeling the protein. The coordinates of heavy atoms from the crystal structure (PDB code 2IYG) of the protein in the dark state served as starting point for 10 ns MD simulations. Representative MD frames were used in QM(B3LYP/cc-pVDZ)/MM(AMBER) calculations to locate minimum energy configurations of the model system. Vertical electronic excitation energies were estimated for the molecular clusters comprising the quantum subsystems of the QM/MM optimized structures using the SOS-CIS(D) quantum chemistry method. Computational results support the occurrence of photoreaction intermediates that are characterized by spectral absorption bands between those of the dark and light states. They agree with crystal structures of reaction intermediates (PDB code 2IYI) observed in the AppA BLUF domain. Transformations of the Gln63 side chain stimulated by photo-excitation and performed with the assistance of the chromophore and the Met106 side chain are responsible for these intermediates.     

Accurate assessment of the strain energy in a protein-bound drug using QM/MM X-ray refinement and converged quantum chemistry

An ongoing question regarding the energetics of protein‐ligand binding has been; what is the strain energy that a ligand pays (if any) when binding to its protein target? The traditional method to estimate strain energy uses force fields to calculate the energy difference between the ligand bound conformation and its nearest local minimum/global minimum on the gas‐phase or aqueous phase potential energy surface. This makes the implicit assumption that the underlying force field as well as the reference crystal structure is accurate. Herein, we use ibuprofen as a test case and compare MMFF and ab initio QM methods to identify the local and global minimum conformations. Nine low energy conformations were identified with HF/6‐31G* geometry optimization in vacuo. We also obtained highly accurate relative energies for ibuprofen's conformational energy surface based on M06/aug‐cc‐pVXZ (X = D and T), MP2/aug‐cc‐pVXZ (X = D and T) and the MP2/CBS method (with and without solvent corrections). Moreover, we curate and re‐refine the ibuprofen‐protein complex (PDB 2BXG) using QM/MM X‐ray refinement approaches (HF/6‐31G* was the QM method and the MM model was the AMBER force field ff99sb), which were compared with the low energy conformers to calculate the strain energy. The result indicates that there was an 88% reduction in ibuprofen conformation strain using the QM/MM refined structure versus the original PDB ibuprofen conformations. Furthermore, our results indicate that, due to its inherent limitations in estimating electrostatic interactions, force fields are not suitable to gauge strain energy for charged drug molecules like ibuprofen. The present work offers a carefully validated conformational potential energy surface for a drug molecule as well as a reliable QM/MM re‐refined X‐ray structure that can be used to test current structure‐based drug design approaches. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

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     

Molecular dynamics and free energy perturbation study of hydride-ion transfer step in dihydrofolate reductase using combined quantum and molecular mechanical model

We used molecular dynamics simulation and free energy perturbation (FEP) methods to investigate the hydride-ion transfer step in the mechanism for the nicotinamide adenine dinucleotide phosphate (NADPH)-dependent reduction of a novel substrate by the enzyme dihydrofolate reductase (DHFR). The system is represented by a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1 semiempirical molecular orbital method for the reacting substrate and NADPH cofactor fragments, the AMBER force field for DHFR, and the TIP3P model for solvent water. The FEP calculations were performed for a number of choices for the QM system. The substrate, 8-methylpterin, was treated quantum mechanically in all the calculations, while the larger cofactor molecule was partitioned into various QM and MM regions with the addition of “link” atoms (F, CH₃, and H). Calculations were also carried out with the entire NADPH molecule treated by QM. The free energies of reaction and the net charges on the NADPH fragments were used to determine the most appropriate QM/MM model. The hydride-ion transfer was also carried out over several FEP pathways, and the QM and QM/MM component free energies thus calculated were found to be state functions (i.e., independent of pathway). A ca. 10 kcal/mol increase in free energy for the hydride-ion transfer with an activation barrier of ca. 30 kcal/mol was calculated. The increase in free energy on the hydride-ion transfer arose largely from the QM/MM component. Analysis of the QM/MM energy components suggests that, although a number of charged residues may contribute to the free energy change through long-range electrostatic interactions, the only interaction that can account for the 10 kcal/mol increase in free energy is the hydrogen bond between the carboxylate side chain of Glu30 (avian DHFR) and the activated (protonated) substrate. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 977–988, 1998