Sampling free energy surfaces as slices by combining umbrella sampling and metadynamics

Metadynamics (MTD) is a very powerful technique to sample high‐dimensional free energy landscapes, and due to its self‐guiding property, the method has been successful in studying complex reactions and conformational changes. MTD sampling is based on filling the free energy basins by biasing potentials and thus for cases with flat, broad, and unbound free energy wells, the computational time to sample them becomes very large. To alleviate this problem, we combine the standard Umbrella Sampling (US) technique with MTD to sample orthogonal collective variables (CVs) in a simultaneous way. Within this scheme, we construct the equilibrium distribution of CVs from biased distributions obtained from independent MTD simulations with umbrella potentials. Reweighting is carried out by a procedure that combines US reweighting and Tiwary–Parrinello MTD reweighting within the Weighted Histogram Analysis Method (WHAM). The approach is ideal for a controlled sampling of a CV in a MTD simulation, making it computationally efficient in sampling flat, broad, and unbound free energy surfaces. This technique also allows for a distributed sampling of a high‐dimensional free energy surface, further increasing the computational efficiency in sampling. We demonstrate the application of this technique in sampling high‐dimensional surface for various chemical reactions using ab initio and QM/MM hybrid molecular dynamics simulations. Further, to carry out MTD bias reweighting for computing forward reaction barriers in ab initio or QM/MM simulations, we propose a computationally affordable approach that does not require recrossing trajectories. © 2016 Wiley Periodicals, Inc.     

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.     

Serine protease acylation proceeds with a subtle re-orientation of the histidine ring at the tetrahedral intermediate

The acylation mechanism of a prototypical serine protease trypsin and its complete free energy reaction profile have been determined by Born-Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling.     

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.     

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.     

Multiscale electrostatic embedding simulations for modeling structure and dynamics of molecules in solution: A tutorial review

The main concepts and important technical details of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) simulations are explained and illustrated with the intent of assisting newcomers in performing and gauging the accuracy of such simulations, focused on smaller molecules in solution. Beginners are advised on how to increase the reliability and accuracy of the simulations through benchmarking. Central considerations on methodologies for QM/MM Molecular Dynamics (MD) simulations are presented, alongside technical fundamentals regarding the construction and manipulation of simulation systems using the python-based Atomic Simulation Environment (ASE). A worked example of QM/MM Born–Oppenheimer MD is included, and a flowchart summarizing the most salient decisions and tasks within the methodology is presented.     

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     

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.     

Ab initio simulations of two‐dimensional electronic spectra: The SOS//QM/MM approach

Two‐dimensional electronic spectroscopy (2DES) is a cutting‐edge technique for investigating with high temporal resolution energy transfer, structure, and dynamics in a wide range of systems in physical chemistry, energy sciences, biophysics, and biocatalysis. However, the interpretation of 2DES is challenging and requires computational modeling. This perspective provides a roadmap for the development of computational tools that could be routinely applied to simulate 2DES spectra of multichromophoric systems active in the UV region (2DUV) using state‐of‐the‐art ab initio electronic structure methods within a quatum mechanics/molecular mechanics (QM/MM) scheme and the sum‐over‐states (SOS) approach (here called SOS//QM/MM). Multiconfigurational and multireference perturbative methods, such as the complete active space self‐consistent field and second‐order multireference perturbation theory (CASPT2) techniques, can be applied to reliably calculate the electronic properties of multichromophoric systems. Hybrid QM/MM method and molecular dynamics techniques can be used to assess environmental and conformational effects, respectively, that shape the 2D electronic spectra. DNA and proteins are important biological targets containing UV chromophores. We report ab initio simulation of 2DUV spectra of a cyclic tetrapeptide containing two interacting aromatic side chains, a model system for the study of protein structure and dynamics by means of 2DUV spectroscopy. © 2013 Wiley Periodicals, Inc.     

Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: sequential sampling and optimization on the potential of mean force surface

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient. In our new sequential sampling and optimization approach, we aim to reduce the amount of MM sampling while still retaining the accuracy of the results by first carrying out MM phase-space sampling and then optimizing the QM subsystem in the fixed-size ensemble of MM conformations. The resulting QM optimized structures are then used to obtain more accurate sampling of the MM subsystem. This process of sequential MM sampling and QM optimization is iterated until convergence. The use of a fixed-size, finite MM conformational ensemble enables the precise evaluation of the QM potential of mean force and its gradient within the ensemble, thus circumventing the challenges associated with statistical averaging and significantly speeding up the convergence of the optimization process. To further improve the accuracy of the QM/MM-MFEP method, the reaction path potential method developed by Lu and Yang [Z. Lu and W. Yang, J. Chem. Phys. 121, 89 (2004)] is employed to describe the QM/MM electrostatic interactions in an approximate yet accurate way with a computational cost that is comparable to classical MM simulations. The new method was successfully applied to two example reaction processes, the classical SN2 reaction of Cl-+CH3Cl in solution and the second proton transfer step of the reaction catalyzed by the enzyme 4-oxalocrotonate tautomerase. The activation free energies calculated with this new sequential sampling and optimization approach to the QM/MM-MFEP method agree well with results from other simulation approaches such as the umbrella sampling technique with direct QM/MM dynamics sampling, demonstrating the accuracy of the iterative QM/MM-MFEP method.     

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.     

Adjustment of Born‐Oppenheimer electronic wave functions to simplify close coupling calculations

Technical problems connected with use of the Born‐Oppenheimer clamped‐nuclei approximation to generate electronic wave functions, potential energy surfaces (PES), and associated properties are discussed. A computational procedure for adjusting the phases of the wave functions, as well as their order when potential crossings occur, is presented which is based on the calculation of overlaps between sets of molecular orbitals and configuration interaction eigenfunctions obtained at neighboring nuclear conformations. This approach has significant advantages for theoretical treatments describing atomic collisions and photo‐dissociation processes by means of ab initio PES, electronic transition moments, and nonadiabatic radial and rotational coupling matrix elements. It ensures that the electronic wave functions are continuous over the entire range of nuclear conformations considered, thereby greatly simplifying the process of obtaining the above quantities from the results of single‐point Born‐Oppenheimer calculations. The overlap results are also used to define a diabatic transformation of the wave functions obtained for conical intersections that greatly simplifies the computation of off‐diagonal matrix elements by eliminating the need for complex phase factors. © 2013 Wiley Periodicals, Inc.     

A molecular mechanics model for imatinib and imatinib:kinase binding

Imatinib is an important anticancer drug, which binds specifically to the Abl kinase and blocks its signalling activity. To model imatinib:protein interactions, we have developed a molecular mechanics force field for imatinib and four close analogues, which is consistent with the CHARMM force field for proteins and nucleic acids. Atomic charges and Lennard‐Jones parameters were derived from a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule interacting with imatinib fragments at 32 different positions. We considered both a neutral and a protonated imatinib. The final RMS deviation between the ab initio and force field energies, averaged over both forms, was 0.2 kcal/mol. The model also reproduces the ab initio geometry and flexibility of imatinib. To apply the force field to imatinib:Abl simulations, it is also necessary to determine the most likely imatinib protonation state when it binds to Abl. This was done using molecular dynamics free energy simulations, where imatinib is reversibly protonated during a series of MD simulations, both in solution and in complex with Abl. The simulations indicate that imatinib binds to Abl in its protonated, positively‐charged form. To help test the force field and the protonation prediction, we did MD free energy simulations that compare the Abl binding affinities of two imatinib analogs, obtaining good agreement with experiment. Finally, two new imatinib variants were considered, one of which is predicted to have improved Abl binding. This variant could be of interest as a potential drug. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010