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.     

Theoretical study of vitamin properties from combined QM-MM methods: Comparison of chemical shifts and energy

The combination of quantum mechanics (QM) and molecular mechanics (MM) methods has become an alternative tool for many applications for which pure QM and MM are not suitable. The QM-MM method has been used for different types of problems, for example, structural biology, surface phenomena, and the liquid phase. In this paper, we have implemented these methods for vitamins, an important kind of biological molecule, and then compared results. The calculations were done by the full ab initio method (HF/3–21 g and HF/6–31 g) and QM-MM (ONIOM) method with HF(3–21 g)/AM1/UFF; then, we found that the geometry obtained by the QM-MM method is very accurate and this rapid method can be used in place of time consuming ab initio methods for large molecules. A comparison of energy values in the QM-MM and QM methods is given. We compare chemical shifts and conclude that the QM-MM method is a perturbed full QM method. The text was submitted by the authors in English.     

MSCALE: A General Utility for Multiscale Modeling

The combination of theoretical models of macromolecules that exist at different spatial and temporal scales has become increasingly important for addressing complex biochemical problems. This work describes the extension of concurrent multiscale approaches, introduces a general framework for carrying out calculations, and describes its implementation into the CHARMM macromolecular modeling package. This functionality, termed MSCALE, generalizes both the additive and subtractive multiscale scheme (e.g. QM/MM ONIOM-type), and extends its support to classical force fields, coarse grained modeling (e.g. ENM, GNM, etc.), and a mixture of them all. The MSCALE scheme is completely parallelized with each subsystem running as an independent, but connected calculation. One of the most attractive features of MSCALE is the relative ease of implementation using the standard MPI communication protocol. This allows external access to the framework and facilitates the combination of functionality previously isolated in separate programs. This new facility is fully integrated with free energy perturbation methods, Hessian based methods, and the use of periodicity and symmetry, which allows the calculation of accurate pressures. We demonstrate the utility of this new technique with four examples; (1) subtractive QM/MM and QM/QM calculations; (2) multi-force field alchemical free energy perturbation; (3) integration with the SANDER module of AMBER and the TINKER package to gain access to potentials not available in CHARMM; and (4) mixed resolution (i.e. coarse grain / all-atom) normal mode analysis. The potential of this new tool is clearly established and in conclusion an interesting mathematical problem is highlighted and future improvements are proposed.     

Combined quantum mechanical/molecular mechanical study on the pentacoordinated ferric and ferrous cytochrome P450cam complexes

The pentacoordinated ferric and ferrous cytochrome P450(cam) complexes have been investigated by combined quantum mechanical/molecular mechanical (QM/MM) calculations in the presence of a protein/solvent environment and by QM calculations on the isolated QM regions with use of density functional theory. The B3LYP functional has been found more reliable than the BLYP and BHLYP functionals for estimating the relative state energies. The B3LYP/CHARMM calculations with an all-electron basis set for iron give high-spin ground states for the title complexes, in agreement with experiment. The comparison of the B3LYP/CHARMM results of the entire protein system with the B3LYP calculations on the naked QM regions shows that the amount of stabilization by the protein environment is largest for the intermediate-spin states, followed by the high-spin states of the complexes. The calculation of M?ssbauer parameters in the presence of the enzyme environment confirms the double occupation of the d(xz) orbital in the quintet spin state of the ferrous complex, consistent with the computed QM/MM energies in the enzyme environment, while the d(x)2(-)(y)2 orbital is doubly occupied in the gas-phase quintet state.     

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.     

Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations

A hybrid quantum mechanical/molecular mechanical (QM/MM) potential energy function with Hartree-Fock, density functional theory (DFT), and post-HF (RIMP2, MP2, CCSD) capability has been implemented in the CHARMM and Q-Chem software packages. In addition, we have modified CHARMM and Q-Chem to take advantage of the newly introduced replica path and the nudged elastic band methods, which are powerful techniques for studying reaction pathways in a highly parallel (i.e., parallel/parallel) fashion, with each pathway point being distributed to a different node of a large cluster. To test our implementation, a series of systems were studied and comparisons were made to both full QM calculations and previous QM/MM studies and experiments. For instance, the differences between HF, DFT, MP2, and CCSD QM/MM calculations of H2O...H2O, H2O...Na+, and H2O...Cl- complexes have been explored. Furthermore, the recently implemented polarizable Drude water model was used to make comparisons to the popular TIP3P and TIP4P water models for doing QM/MM calculations. We have also computed the energetic profile of the chorismate mutase catalyzed Claisen rearrangement at various QM/MM levels of theory and have compared the results with previous studies. Our best estimate for the activation energy is 8.20 kcal/mol and for the reaction energy is -23.1 kcal/mol, both calculated at the MP2/6-31+G(d)//MP2/6-31+G(d)/C22 level of theory.     

A multicentered approach to integrated QM/QM calculations. Applications to multiply hydrogen bonded systems

A multicentered integrated QM/QM technique has been developed. By separating high-level calculations in distinct regions of molecules, the multicentered approach supplants a single large high-level calculation with several smaller calculations. Due to the steep polynomial scaling of traditional ab initio quantum chemical methods, this separation significantly enhances the computational efficiency of QM/QM methods. The straightforward implementation of this multicentered approach is illustrated with several large poly-alcohols that form hydrogen bonds with water. The largest alcohol-water complex contains 81 atoms. For properly selected model systems, this multicentered approach introduces essentially no error in the dissociation energies of these complexes relative to conventional QM/QM schemes. This multicentered technique should be easily extended to other, more general integrated methods (QM/MM, ONIOM, etc).     

Geometry optimization using tuned and balanced redistributed charge schemes for combined quantum mechanical and molecular mechanical calculations

We performed geometry optimizations using the tuned and balanced redistributed charge algorithms to treat the QM-MM boundary in combined quantum mechanical and molecular mechanical (QM/MM) methods. In the tuned and balanced redistributed charge (TBRC) scheme, the QM boundary atom is terminated by a tuned F link atom, and the charge of the MM boundary atom is properly adjusted to conserve the total charge of the entire QM/MM system; then the adjusted MM boundary charge is moved evenly to the midpoints of the bonds between the MM boundary atom and its neighboring MM atoms. In the tuned and balanced redistributed charge-2 (TBRC2) scheme, the adjusted MM boundary charge is moved evenly to all MM atoms that are attached to the MM boundary atom. A new option, namely charge smearing, has been added to the TBRC scheme, yielding the tuned and balanced smeared redistributed charge (TBSRC) scheme. In the new scheme, the redistributed charges near the QM-MM boundary are smeared to make the electrostatic interactions between the QM region and the redistributed charges more realistic. The TBRC2 scheme and new TBSRC scheme have been tested for various kinds of bonds at a QM-MM boundary, including C-C, C-N, C-O, O-C, N-C, C-S, S-S, S-C, C-Si, and O-N bonds. Charge smearing is necessary if the redistributed charges are close to the QM region, as in the TBSRC scheme, but not if the redistributed charge is farther from the QM region, as in the TBRC2 scheme. We found that QM/MM results using either the TBRC2 scheme or the TBSRC scheme agree well with full QM results; the mean unsigned error (MUE) of the QM/MM deprotonation energy is 1.6 kcal/mol in both cases, and the MUE of QM/MM optimized bond lengths over the three bonds closest to the QM-MM boundary, with errors averaged over the protonated forms and unprotonated forms, is 0.015 ? for TBRC2 and 0.021 ? for TBSRC. The improvements in the new scheme are essential for QM-MM boundaries that pass through a polar bond, but even for boundaries that pass through C-C bonds, the improvement can be quite significant.     

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.     

thermocalc — A poor man's approach to computational thermochemistry

We present thermocalc, a Perl module to perform the automated calculation of atomization energies and heats of formation for lists of molecules. The methods used are based on density functional theory and second‐order perturbation theory to ensure that data sets of medium sized to large molecules can be run at reasonable throughput rates. The quantum chemical calculations are performed using the program package TURBOMOLE in a three‐step protocol. In a first step, a pre‐optimization of the structure and a zero‐point energy calculation are performed. As second step, a geometry optimization is being carried out, and the last step is a single point energy calculation. The level of theory used in the different steps can be modified by the user to allow for customized protocols. The performance of example protocols is investigated on different test sets of molecules. In the course of this work, a simple, but efficient one‐parameter correction term based on the shared electron numbers has been developed, which reduces the error of calculated heats of formation significantly. © 2012 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.     

Comparison of methods for finding saddle points without knowledge of the final states

Within the harmonic approximation to transition state theory, the biggest challenge involved in finding the mechanism or rate of transitions is the location of the relevant saddle points on the multidimensional potential energy surface. The saddle point search is particularly challenging when the final state of the transition is not specified. In this article we report on a comparison of several methods for locating saddle points under these conditions and compare, in particular, the well-established rational function optimization (RFO) methods using either exact or approximate Hessians with the more recently proposed minimum mode following methods where only the minimum eigenvalue mode is found, either by the dimer or the Lanczos method. A test problem involving transitions in a seven-atom Pt island on a Pt(111) surface using a simple Morse pairwise potential function is used and the number of degrees of freedom varied by varying the number of movable atoms. In the full system, 175 atoms can move so 525 degrees of freedom need to be optimized to find the saddle points. For testing purposes, we have also restricted the number of movable atoms to 7 and 1. Our results indicate that if attempting to make a map of all relevant saddle points for a large system (as would be necessary when simulating the long time scale evolution of a thermal system) the minimum mode following methods are preferred. The minimum mode following methods are also more efficient when searching for the lowest saddle points in a large system, and if the force can be obtained cheaply. However, if only the lowest saddle points are sought and the calculation of the force is expensive but a good approximation for the Hessian at the starting position of the search can be obtained at low cost, then the RFO approaches employing an approximate Hessian represent the preferred choice. For small and medium sized systems where the force is expensive to calculate, the RFO approaches employing an approximate Hessian is also the more efficient, but when the force and Hessian can be obtained cheaply and only the lowest saddle points are sought the RFO approach using an exact Hessian is the better choice. These conclusions have been reached based on a comparison of the total computational effort needed to find the saddle points and the number of saddle points found for each of the methods. The RFO methods do not perform very well with respect to the latter aspect, but starting the searches further away from the initial minimum or using the hybrid RFO version presented here improves this behavior considerably in most cases.     

Solvent effects on peroxynitrite structure and properties from QM/MM simulations

We have investigated the structure and the vibrational spectrum of peroxynitrite anion in aqueous solution by means of combined quantum-classical (QM/MM) molecular dynamics simulations. In our QM/MM scheme, the reactant was modeled using density functional theory with a Gaussian basis set and the solvent was described using the mean-field TIP4P and the polarizable TIP4P-FQ force fields. The choice of basis sets, functionals and force field parameters has been validated by performing calculations on isolated peroxynitrite and on small peroxynitrite-water complexes. Poor values for isolated peroxynitrite structural properties and vibrational frequencies are found for most ab initio methods, particularly regarding the ON-OO(-) bond distance and the harmonic stretching frequency, probably due to the singlet-triplet instability found in the HF wave function. On the other hand, DFT methods yield results in better agreement with high level CCSD(T) ab initio calculations. We have computed the vibrational spectrum for aqueous peroxynitrite by calculating the Fourier transform of the velocity autocorrelation function extracted from the QM-MM molecular dynamics simulations. Our computational scheme, which allows for the inclusion of both anharmonicity and solvent effects, is able to clarify previous discrepancies regarding the experimental spectra assignments and to shed light on the subtle interplay between solvation and peroxynitrite structure and properties.     

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     

pKa calculations in solution and proteins with QM/MM free energy perturbation simulations: a quantitative test of QM/MM protocols

The accuracy of biological simulations depends, in large part, on the treatment of electrostatics. Due to the availability of accurate experimental values, calculation of pKa provides stringent evaluation of computational methods. The generalized solvent boundary potential (GSBP) and Ewald summation electrostatic treatments were recently implemented for combined quantum mechanical and molecular mechanics (QM/MM) simulations by our group. These approaches were tested by calculating pKa shifts due to differences in electronic structure and electrostatic environment; the shifts were determined for a series of small molecules in solution, using various electrostatic treatments, and two residues (His 31, Lys 102) in the M102K T4-lysozyme mutant with large pKa shifts, using the GSBP approach. The calculations utilized a free energy perturbation scheme with the QM/MM potential function involving the self-consistent charge density functional tight binding (SCC-DFTB) and CHARMM as the QM and MM methods, respectively. The study of small molecules demonstrated that inconsistent electrostatic models produced results that were difficult to correct in a robust manner; by contrast, extended electrostatics, GSBP, and Ewald simulations produced consistent results once a bulk solvation contribution was carefully chosen. In addition to the electrostatic treatment, the pKa shifts were also sensitive to the level of the QM method and the scheme of treating QM/MM Coulombic interactions; however, simple perturbative corrections based on SCC-DFTB/CHARMM trajectories and higher level single point energy calculations were found to give satisfactory results. Combining all factors gave a root-mean-square difference of 0.7 pKa units for the relative pKa values of the small molecules compared to experiment. For the residues in the lysozyme, an accurate pKa shift was obtained for His 31 with multiple nanosecond simulations. For Lys 102, however, the pKa shift was estimated to be too large, even after more than 10 nanosecond simulations for each lambda window; the difficulty was due to the significant, but slow, reorganization of the protein and water structure when Lys 102 was protonated. The simulations support that Lys 102 is deprotonated in the X-ray structure and the protein is highly destabilized when this residue is protonated.     

Transition structure selectivity in enzyme catalysis: a QM/MM study of chorismate mutase

Two different transition structures (TSs) have been located and characterized for the chorismate conversion to prephenate in Bacillus subtilis chorismate mutase by means of hybrid quantum-mechanical/molecular-mechanical (QM/MM) calculations. GRACE software, combined with an AM1/CHARMM24/TIP3P potential, has been used involving full gradient relaxation of the position of ca. 3300 atoms. These TSs have been connected with their respective reactants and products by the intrinsic reaction coordinate (IRC) procedure carried out in the presence of the protein environment, thus obtaining for the first time a realistic enzymatic reaction path for this reaction. Similar QM/MM computational schemes have been applied to study the chemical reaction solvated by ca. 500 water molecules. Comparison of these results together with gas phase calculations has allowed understanding of the catalytic efficiency of the protein. The enzyme stabilizes one of the TSs (TS_OHout) by means of specific hydrogen bond interactions, while the other TS (TS_OHin) is the preferred one in vacuum and in water. The enzyme TS is effectively more polarized but less dissociative than the corresponding solvent and gas phase TSs. Electrostatic stabilization and an intramolecular charge-transfer process can explain this enzymatically induced change. Our theoretical results provide new information on an important enzymatic transformation and the key factors responsible for efficient selectivity are clarified. Received: 25 March 2000 / Accepted: 7 August 2000 / Published online: 23 November 2000     

Reaction path potential for complex systems derived from combined ab initio quantum mechanical and molecular mechanical calculations

Combined ab initio quantum mechanical and molecular mechanical calculations have been widely used for modeling chemical reactions in complex systems such as enzymes, with most applications being based on the determination of a minimum energy path connecting the reactant through the transition state to the product in the enzyme environment. However, statistical mechanics sampling and reaction dynamics calculations with a combined ab initio quantum mechanical (QM) and molecular mechanical (MM) potential are still not feasible because of the computational costs associated mainly with the ab initio quantum mechanical calculations for the QM subsystem. To address this issue, a reaction path potential energy surface is developed here for statistical mechanics and dynamics simulation of chemical reactions in enzymes and other complex systems. The reaction path potential follows the ideas from the reaction path Hamiltonian of Miller, Handy and Adams for gas phase chemical reactions but is designed specifically for large systems that are described with combined ab initio quantum mechanical and molecular mechanical methods. The reaction path potential is an analytical energy expression of the combined quantum mechanical and molecular mechanical potential energy along the minimum energy path. An expansion around the minimum energy path is made in both the nuclear and the electronic degrees of freedom for the QM subsystem internal energy, while the energy of the subsystem described with MM remains unchanged from that in the combined quantum mechanical and molecular mechanical expression and the electrostatic interaction between the QM and MM subsystems is described as the interaction of the MM charges with the QM charges. The QM charges are polarizable in response to the changes in both the MM and the QM degrees of freedom through a new response kernel developed in the present work. The input data for constructing the reaction path potential are energies, vibrational frequencies, and electron density response properties of the QM subsystem along the minimum energy path, all of which can be obtained from the combined quantum mechanical and molecular mechanical calculations. Once constructed, it costs much less for its evaluation. Thus, the reaction path potential provides a potential energy surface for rigorous statistical mechanics and reaction dynamics calculations of complex systems. As an example, the method is applied to the statistical mechanical calculations for the potential of mean force of the chemical reaction in triosephosphate isomerase.