Calibration of the Quantum/Classical Hamiltonian in Semiempirical QM/MM AM1 and PM3 Methods

Minimum Free Energy Pathways of Reactive Processes with Nudged Elastic Bands

The determination of minimum free energy pathways (MFEP) is one of the most widely used strategies to study reactive processes. For chemical reactions in complex environments, the combination of quantum mechanics (QM) with a molecular mechanics (MM) representation is usually necessary in a hybrid QM/MM framework. However, even within the QM/MM approximation, the affordable sampling of the phase space is, in general, quite restricted. To reduce drastically the computational cost of the simulations, several methods such as umbrella sampling require performing a priori a selection of a reaction coordinate. The quality of the computed results, in an affordable computational time, is intimately related to the reaction coordinate election which is, in general, a nontrivial task. In this work, we provide an approach to model reactive processes in complex environments that does not require the a priori selection of a reaction coordinate. The proposed methodology combines QM/MM simulations with an extrapolation of the nudged elastic bands (NEB) method to the free energy surface (FENEB). We present and apply our own FENEB scheme to optimize MFEP in different reactive processes, using QM/MM frameworks at semiempirical and density functional theory levels. Our implementation is based on performing the FENEB optimization by uncoupling the optimization of the band in a perpendicular and tangential direction. In each step, a full optimization with the spring force is performed, which guarantees that the images remain evenly distributed. The robustness of the method and the influence of sampling on the quality of the optimized MFEP and its associated free energy barrier are studied. We show that the FENEB method provides a good estimation of the reaction barrier even with relatively short simulation times, supporting that its combination with QM/MM frameworks provides an adequate tool to study chemical processes in complex environments.

Semi-Empirical Born-Oppenheimer Molecular Dynamics (SEBOMD) within the Amber Biomolecular Package

Semi-empirical quantum methods from the neglect of differential diatomic overlap (NDDO) family such as MNDO, AM1, or PM3 are fast albeit approximate quantum methods. By combining them with linear scaling methods like the divide & conquer (D&C) method, it is possible to quickly evaluate the energy of systems containing hundreds to thousands of atoms. We here present our implementation in the Amber biomolecular package of a SEBOMD module that provides a way to run semi-empirical Born-Oppenheimer molecular dynamics. At each step of a SEBOMD, a fully converged self-consistent field (SCF) calculation is performed to obtain the semiempirical quantum potential energy of a molecular system encaged or not in periodic boundary conditions. We describe the implementation and the features of our SEBOMD implementation. We show the requirements to conserve the total energy in NVE simulations, and how to accelerate SCF convergence through density matrix extrapolation. Specific ways of handling periodic boundary conditions using mechanical embedding or electrostatic embedding through a tailored quantum Ewald summation is developed. The parallel performance of SEBOMD simulations using the D&C scheme are presented for liquid water systems of various sizes, and a comparison between the traditional full diagonalization scheme and the D&C approach for the reproduction of the structure of liquid water illustrates the potentiality of SEBOMD to simulate molecular systems containing several hundreds of atoms for hundreds of picoseconds with a quantum mechanical potential in a reasonable amount of CPU time.

A combination of ternary classification models and reporter gene assays for the comprehensive thyroid hormone disruption profiles of 209 polychlorinated biphenyls

Computational toxicology is widely applied to screen tens and thousands of potential environmental endocrine disruptors (EDCs) for its great efficiency and rapid evaluation in recent years. Polychlorinated biphenyls (PCBs) with 209 congeners have not been comprehensively tested for their ability to interact with the thyroid receptor (TR), which is one of the most extensively studied targets related to the effects of EDCs. In this study, we aimed to determine the thyroid-disrupting mechanisms of 209 PCBs through the combination of a novel computational ternary classification model and luciferase reporter gene assay. The reporter gene assay was performed to examine the hormone activities of 22 PCBs via TR to classify their profiles as agonistic, antagonistic or inactive. Thus, six agonists, eleven antagonists and seven inactive chemicals against TR were identified in in vitro assays. Then, six relevant variables, including molecular structural descriptors and molecular docking scores, were selected for describing compounds. Machine learning methods (i.e., linear discriminant analysis (LDA) and support vector machines (SVM)) were used to build classification models. The optimal model was obtained by using SVM, with an accuracy of 88.24% in the training set, 80.0% in the test set and 75.0% in 10-fold cross-validation. The remaining 187 PCB congeners' hormone activities were predicted using the obtained models. Out of these PCBs, the SVM model predicted 38 agonists and 81 antagonists. The findings revealed that higher chlorinated PCBs tended to have TR-antagonistic activities, whereas lower chlorinated PCBs were agonists. This study provided a reference for further exploring PCBs' thyroid effect.

The influence of the metal cations and microhydration on the reaction trajectory of the N3 ↔ O2 thymine proton transfer: Quantum mechanical study

This study involves the intramolecular proton transfer (PT) process on a thymine nucleobase between N3 and O2 atoms. We explore a mechanism for the PT assisted by hexacoordinated divalent metals cations, namely Mg²⁺ , Zn²⁺ , and Hg²⁺ . Our results point out that this reaction corresponds to a two-stage process. The first involves the PT from one of the aqua ligands toward O2. The implications of this stage are the formation of a hydroxo anion bound to the metal center and a positively charged thymine. To proceed to the second stage, a structural change is needed to allow the negatively charged hydroxo ligand to abstract the N3 proton, which represents the final product of the PT reaction. In the presence of the selected hexaaqua cations, the activation barrier is at most 8 kcal/mol. © 2017 Wiley Periodicals, Inc.

Reaching multi-nanosecond timescales in combined QM/MM molecular dynamics simulations through parallel horsetail sampling

We report an enhanced sampling technique that allows to reach the multi-nanosecond timescale in quantum mechanics/molecular mechanics molecular dynamics simulations. The proposed technique, called horsetail sampling, is a specific type of multiple molecular dynamics approach exhibiting high parallel efficiency. It couples a main simulation with a large number of shorter trajectories launched on independent processors at periodic time intervals. The technique is applied to study hydrogen peroxide at the water liquid-vapor interface, a system of considerable atmospheric relevance. A total simulation time of a little more than 6 ns has been attained for a total CPU time of 5.1 years representing only about 20 days of wall-clock time. The discussion of the results highlights the strong influence of the solvation effects at the interface on the structure and the electronic properties of the solute. © 2017 Wiley Periodicals, Inc.

A polarizable QM/MM approach to the molecular dynamics of amide groups solvated in water

The infrared (IR) spectra of polypeptides are dominated by the so-called amide bands. Because they originate from the strongly polar and polarizable amide groups (AGs) making up the backbone, their spectral positions sensitively depend on the local electric fields. Aiming at accurate computations of these IR spectra by molecular dynamics (MD) simulations, which derive atomic forces from a hybrid quantum and molecular mechanics (QM/MM) Hamiltonian, here we consider the effects of solvation in bulk liquid water on the amide bands of the AG model compound N-methyl-acetamide (NMA). As QM approach to NMA we choose grid-based density functional theory (DFT). For the surrounding MM water, we develop, largely based on computations, a polarizable molecular mechanics (PMM) model potential called GP6P, which features six Gaussian electrostatic sources (one induced dipole, five static partial charge distributions) and, therefore, avoids spurious distortions of the DFT electron density in hybrid DFT/PMM simulations. Bulk liquid GP6P is shown to have favorable properties at the thermodynamic conditions of the parameterization and beyond. Lennard-Jones (LJ) parameters of the DFT fragment NMA are optimized by comparing radial distribution functions in the surrounding GP6P liquid with reference data obtained from a "first-principles" DFT-MD simulation. Finally, IR spectra of NMA in GP6P water are calculated from extended DFT/PMM-MD trajectories, in which the NMA is treated by three different DFT functionals (BP, BLYP, B3LYP). Method-specific frequency scaling factors are derived from DFT-MD simulations of isolated NMA. The DFT/PMM-MD simulations with GP6P and with the optimized LJ parameters then excellently predict the effects of aqueous solvation and deuteration observed in the IR spectra of NMA. As a result, the methods required to accurately compute such spectra by DFT/PMM-MD also for larger peptides in aqueous solution are now at hand.

Application of a BOSS-Gaussian interface for QM/MM simulations of Henry and methyl transfer reactions

Hybrid quantum mechanics and molecular mechanics (QM/MM) computer simulations have become an indispensable tool for studying chemical and biological phenomena for systems too large to treat with QM alone. For several decades, semiempirical QM methods have been used in QM/MM simulations. However, with increased computational resources, the introduction of ab initio and density function methods into on-the-fly QM/MM simulations is being increasingly preferred. This adaptation can be accomplished with a program interface that tethers independent QM and MM software packages. This report introduces such an interface for the BOSS and Gaussian programs, featuring modification of BOSS to request QM energies and partial atomic charges from Gaussian. A customizable C-shell linker script facilitates the interprogram communication. The BOSS-Gaussian interface also provides convenient access to Charge Model 5 (CM5) partial atomic charges for multiple purposes including QM/MM studies of reactions. In this report, the BOSS-Gaussian interface is applied to a nitroaldol (Henry) reaction and two methyl transfer reactions in aqueous solution. Improved agreement with experiment is found by determining free-energy surfaces with MP2/CM5 QM/MM simulations than previously reported investigations using semiempirical methods.

Water interactions with hydrophobic groups: assessment and recalibration of semiempirical molecular orbital methods

In this work, we present a study of the ability of different semiempirical methods to describe intermolecular interactions in water solution. In particular, we focus on methods based on the Neglect of Diatomic Differential Overlap approximation. Significant improvements of these methods have been reported in the literature in the past years regarding the description of non-covalent interactions. In particular, a broad range of methodologies has been developed to deal with the properties of hydrogen-bonded systems, with varying degrees of success. In contrast, the interactions between water and a molecule containing hydrophobic groups have been little analyzed. Indeed, by considering the potential energy surfaces obtained using different semiempirical Hamiltonians for the intermolecular interactions of model systems, we found that none of the available methods provides an entirely satisfactory description of both hydrophobic and hydrophilic interactions in water. In addition, a vibrational analysis carried out in a model system for these interactions, a methane clathrate cluster, showed that some recent methods cannot be used to carry out studies of vibrational properties. Following a procedure established in our group [M. I. Bernal-Uruchurtu, M. T. C. Martins-Costa, C. Millot, and M. F. Ruiz-López, J. Comput. Chem. 21, 572 (2000); W. Harb, M. I. Bernal-Uruchurtu, and M. F. Ruiz-López, Theor. Chem. Acc. 112, 204 (2004)], we developed new parameters for the core-core interaction terms based on fitting potential energy curves obtained at the MP2 level for our model system. We investigated the transferability of the new parameters to describe a system, having both hydrophilic and hydrophobic groups, interacting with water. We found that only by introducing two different sets of parameters for hydrophilic and hydrophobic hydrogen atom types we are able to match the features of the ab initio calculated properties. Once this assumption is made, a good agreement with the MP2 reference is achieved. The results reported in this work provide therefore a direction for future developments of semiempirical approaches that are still required to investigate chemical processes in biomolecules and in large disordered systems.

Multiscale Free Energy Simulations: An Efficient Method for Connecting Classical MD Simulations to QM or QM/MM Free Energies Using Non-Boltzmann Bennett Reweighting Schemes

The reliability of free energy simulations (FES) is limited by two factors: (a) the need for correct sampling and (b) the accuracy of the computational method employed. Classical methods (e.g., force fields) are typically used for FES and present a myriad of challenges, with parametrization being a principle one. On the other hand, parameter-free quantum mechanical (QM) methods tend to be too computationally expensive for adequate sampling. One widely used approach is a combination of methods, where the free energy difference between the two end states is computed by, e.g., molecular mechanics (MM), and the end states are corrected by more accurate methods, such as QM or hybrid QM/MM techniques. Here we report two new approaches that significantly improve the aforementioned scheme; with a focus on how to compute corrections between, e.g., the MM and the more accurate QM calculations. First, a molecular dynamics trajectory that properly samples relevant conformational degrees of freedom is generated. Next, potential energies of each trajectory frame are generated with a QM or QM/MM Hamiltonian. Free energy differences are then calculated based on the QM or QM/MM energies using either a non-Boltzmann Bennett approach (QM-NBB) or non-Boltzmann free energy perturbation (NB-FEP). Both approaches are applied to calculate relative and absolute solvation free energies in explicit and implicit solvent environments. Solvation free energy differences (relative and absolute) between ethane and methanol in explicit solvent are used as the initial test case for QM-NBB. Next, implicit solvent methods are employed in conjunction with both QM-NBB and NB-FEP to compute absolute solvation free energies for 21 compounds. These compounds range from small molecules such as ethane and methanol to fairly large, flexible solutes, such as triacetyl glycerol. Several technical aspects were investigated. Ultimately some best practices are suggested for improving methods that seek to connect MM to QM (or QM/MM) levels of theory in FES.

Effects of Dispersion in Density Functional Based Quantum Mechanical/Molecular Mechanical Calculations on Cytochrome P450 Catalyzed Reactions

Density functional theory (DFT) based quantum mechanical/molecular mechanical (QM/MM) calculations have provided valuable insight into the reactivity of the cytochrome P450 family of enzymes (P450s). A failure of commonly used DFT methods, such as B3LYP, is the neglect of dispersion interactions. An empirical dispersion correction has been shown to improve the accuracy of gas phase DFT calculations of P450s. The current work examines the effect of the dispersion correction in QM/MM calculations on P450s. The hydrogen abstraction from camphor, and hydrogen abstraction and C-O addition of cyclohexene and propene by P450cam have been modeled, along with the addition of benzene to Compound I in CYP2C9, at the B3LYP-D2/CHARMM27 level of theory. Single point energy calculations were also performed at the B3LYP-D3//B3LYP-D2/CHARMM27 level. The dispersion corrections lower activation energy barriers significantly (by ∼5 kcal/mol), as seen for gas phase calculations, but has a small effect on optimized geometries.These effects are likely to be important in modeling reactions catalyzed by other enzymes also. Given the low computational cost of including such dispersion corrections, we recommend doing so in all B3LYP based QM/MM calculations.

Analysis of polarization in QM/MM modelling of biologically relevant hydrogen bonds

Combined quantum mechanics/molecular mechanics (QM/MM) methods are increasingly important for the study of chemical reactions and systems in condensed phases. Here, we have tested the accuracy of a density functional theory-based QM/MM implementation (B3LYP/6-311+G(d,p)/CHARMM27) on a set of biologically relevant interactions by comparison with full QM calculations. Intermolecular charge transfer due to hydrogen bond formation is studied to assess the severity of spurious polarization of QM atoms by MM point charges close to the QM/MM boundary. The changes in total electron density and natural bond orbital atomic charges due to hydrogen bond formation in selected complexes obtained at the QM/MM level are compared with full QM results. It is found that charge leakage from the QM atoms to MM atomic point charges close to the QM/MM boundary is not a serious problem, at least with limited basis sets. The results are encouraging in showing that important properties of key biomolecular interactions can be treated well at the QM/MM level employing good-quality levels of QM theory.

Extension of the MST continuum solvation model to the RM1 semiempirical Hamiltonian

The need to simulate large-sized molecules or to deal with large series of compounds is a challenging topic in computational chemistry, which has stimulated the development of accurate semiempirical methods, such as the recently reported Recife Model 1 (RM1; J Comput Chem 2006, 27, 1101). Even though RM1 may prove to be of value simply due to the enhanced quantitative accuracy in gas phase, it is unclear how the new parameters optimized for RM1 affect the suitability of this semiempirical Hamiltonian to study chemical processes in condensed phases. To address this question, we report the parametrization of the MST/RM1 continuum model for neutral solutes in water, octanol, chloroform and carbon tetrachloride, and for ions in water. The results are used to discuss the transferability of the solvation parameters implemented in previous MST/AM1 and MST/PM3 models.

Charge-dependent model for many-body polarization, exchange, and dispersion interactions in hybrid quantum mechanical/molecular mechanical calculations

This work explores a new charge-dependent energy model consisting of van der Waals and polarization interactions between the quantum mechanical (QM) and molecular mechanical (MM) regions in a combined QMMM calculation. van der Waals interactions are commonly treated using empirical Lennard-Jones potentials, whose parameters are often chosen based on the QM atom type (e.g., based on hybridization or specific covalent bonding environment). This strategy for determination of QMMM nonbonding interactions becomes tedious to parametrize and lacks robust transferability. Problems occur in the study of chemical reactions where the "atom type" is a complex function of the reaction coordinate. This is particularly problematic for reactions, where atoms or localized functional groups undergo changes in charge state and hybridization. In the present work we propose a new model for nonelectrostatic nonbonded interactions in QMMM calculations that overcomes many of these problems. The model is based on a scaled overlap model for repulsive exchange and attractive dispersion interactions that is a function of atomic charge. The model is chemically significant since it properly correlates atomic size, softness, polarizability, and dispersion terms with minimal one-body parameters that are functions of the atomic charge. Tests of the model are examined for rare-gas interactions with neutral and charged atoms in order to demonstrate improved transferability. The present work provides a new framework for modeling QMMM interactions with improved accuracy and transferability.

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.

Transition-State Characterization of the Ammonia Ionization Process in Aqueous Solution via the Free-Energy Gradient Method

For the ionization process of ammonia in aqueous solution, the transition-state (TS) structure was fully optimized for the first time on the free-energy surface (FES) by applying the free-energy gradient (FEG) method combined with a hybrid quantum mechanical and molecular mechanical molecular dynamics (QM/MM-MD) method. In aqueous solution, the ionization process was found to proceed by way of a clear TS (R(N1-H5) = 1.512 A), which does not exist in the gas phase. The free-energy (FE) of activation for ionization obtained was 14.7 kcal/mol, within the classical approximation, via the QM/MM-MD FEG method, and is found to be in good agreement with 9.57 kcal/mol estimated from the TS theory using the experimental value of the rate constant. Apart from the dynamic correction, it is indicated that the theoretical value would be improved to be 10.28 kcal/mol if the electronic-state calculation could be executed at the B3LYP/6-31G(d) level of theory.

Are AM1 ligand-protein binding enthalpies good enough for use in the rational design of new drugs?

We have examined the performance of semiempirical quantum mechanical methods in solving the problem of accurately predicting protein-ligand binding energies and geometries. Firstly, AM1 and PM3 geometries and binding enthalpies between small molecules that simulate typical ligand-protein interactions were compared with high level quantum mechanical techniques that include electronic correlation (e.g., MP2 or B3LYP). Species studied include alkanes, aromatic systems, molecules including groups with hypervalent sulfur or with donor or acceptor hydrogen bonding capability, as well as ammonium or carboxylate ions. B3LYP/6-311+G(2d,p) binding energies correlated very well with the BSSE corrected MP2/6-31G(d) values. AM1 binding enthalpies also showed good correlation with MP2 values, and their systematic deviation is acceptable when enthalpies are used for the comparison of interaction energies between ligands and a target. PM3 otherwise gave erratic energy differences in comparison to the B3LYP or MP2 approaches. As one would expect, the geometries of the binding complexes showed the known limitations of the semiempirical and DFT methods. AM1 calculations were subsequently applied to a test set consisting of "real" protein active site-ligand complexes. Preliminary results indicate that AM1 could be a valuable tool for the design of new drugs using proteins as templates. This approach also has a reasonable computational cost. The ligand-protein X-ray structures were reasonably reproduced by AM1 calculations and the corresponding AM1 binding enthalpies are in agreement with the results from the "small molecules" test set.

Lennard-Jones parameters for the combined QM/MM method using the B3LYP/6-31G*/AMBER potential

A combined DFT quantum mechanical and AMBER molecular mechanical potential (QM/MM) is presented for use in molecular modeling and molecular simulations of large biological systems. In our approach we evaluate Lennard-Jones parameters describing the interaction between the quantum mechanical (QM) part of a system, which is described at the B3LYP/6-31+G* level of theory, and the molecular mechanical (MM) part of the system, described by the AMBER force field. The Lennard-Jones parameters for this potential are obtained by calculating hydrogen bond energies and hydrogen bond geometries for a large set of bimolecular systems, in which one hydrogen bond monomer is described quantum mechanically and the other is treated molecular mechanically. We have investigated more than 100 different bimolecular systems, finding very good agreement between hydrogen bond energies and geometries obtained from the combined QM/MM calculations and results obtained at the QM level of theory, especially with respect to geometry. Therefore, based on the Lennard-Jones parameters obtained in our study, we anticipate that the B3LYP/6-31+G*/AMBER potential will be a precise tool to explore intermolecular interactions inside a protein environment.

Structures and stability of salt-bridge in aqueous solution

Structures and stability of salt-bridges in aqueous solutions were investigated using a complex formed from the guanidinium (Gdm+) and formate (FmO-) ions as a model system. The Test-particle model (T-model) potentials to describe the interactions in the Gdm+-H2O, FmO(-)-H2O and Gdm+-FmO- complexes were constructed, tested and applied in molecular dynamics (MD) simulations of the aqueous solutions at 298 K. The three-dimensional structures and energetic of the hydrogen bond (H-bond) networks of water in the first hydration shells of the Gdm+ and FmO- ions, as well as the Gdm+-FmO- complex, were visualized and analyzed using various probability distribution (PD) maps. The structures of the average potential energy landscapes at the H-bond networks were employed to characterize the stability and dynamic behavior of water molecules in the first hydration shells of the solutes. It was observed that water molecules in the first hydration shell of the close-contact Gdm+-FmO- complex form associated H-bond networks, which introduce a net stabilization effect to the ion-pair, whereas those in the interstitial H-bond network destabilize and break the solvent-separated Gdm+-FmO- complex. The present results showed that, in order to provide complete insights into the structures and stability of ion-pairs in aqueous solutions, explicit water molecules have to be included in the model calculations.

Benchmark calculations of proton affinities and gas-phase basicities of molecules important in the study of biological phosphoryl transfer

Benchmark calculations of proton affinities and gas-phase basicities of molecules most relevant to biological phosphoryl transfer reactions are presented and compared with available experimental results. The accuracy of proton affinity and gas-phase basicity results obtained from several multi-level model chemistries (CBS-QB3, G3B3, and G3MP2B3) and density-functional quantum models (PBE0, B1B95, and B3LYP) are assessed and compared. From these data, a set of empirical bond enthalpy, entropy, and free energy corrections are introduced that considerably improve the accuracy and predictive capability of the methods. These corrections are applied to the prediction of proton affinity and gas-phase basicity values of important biological phosphates and phosphoranes for which experimental data does not currently exist. Comparison is made with results from semiempirical quantum models that are commonly employed in hybrid quantum mechanical/molecular mechanical simulations. Data suggest that the design of improved semiempirical quantum models with increased accuracy for relative proton affinity values is necessary to obtain quantitative accuracy for phosphoryl transfer reactions in solution, enzymes, and ribozymes.

Computational methods for the study of enzymic reaction mechanisms III: A perturbation plus QM/MM approach for calculating relative free energies of protonation

We describe a coupling parameter, that is, perturbation, approach to effectively create and annihilate atoms in the quantum mechanical Hamiltonian within the closed shell restricted Hartree-Fock formalism. This perturbed quantum mechanical atom (PQA) method is combined with molecular mechanics (MM) methods (PQA/MM) within a molecular dynamics simulation, to model the protein environment (MM region) effects that also make a contribution to the overall free energy change. Using the semiempirical PM3 method to model the QM region, the application of this PQA/MM method is illustrated by calculation of the relative protonation free energy of the conserved OD2 (Asp27) and the N5 (dihydrofolate) proton acceptor sites in the active site of Escherichia coli dihydrofolate reductase (DHFR) with the bound nicotinamide adenine dinucleotide phosphate (NADPH) cofactor. For a number of choices for the QM region, the relative protonation free energy was calculated as the sum of contributions from the QM region and the interaction between the QM and MM regions via the thermodynamic integration (TI) method. The results demonstrate the importance of including the whole substrate molecule in the QM region, and the overall protein (MM) environment in determining the relative stabilities of protonation sites in the enzyme active site. The PQA/MM free energies obtained by TI were also compared with those estimated by a less computationally demanding nonperturbative method based on the linear response approximation (LRA). For some choices of QM region, the total free energies calculated using the LRA method were in very close agreement with the PQA/MM values. However, the QM and QM/MM component free energies were found to differ significantly between the two methods.

Extension of the PDDG/PM3 and PDDG/MNDO semiempirical molecular orbital methods to the halogens

The new semiempirical methods, PDDG/PM3 and PDDG/MNDO, have been parameterized for halogens. For comparison, the original MNDO and PM3 were also reoptimized for the halogens using the same training set; these modified methods are referred to as MNDO' and PM3'. For 442 halogen-containing molecules, the smallest mean absolute error (MAE) in heats of formation is obtained with PDDG/PM3 (5.6 kcal/mol), followed by PM3' (6.1 kcal/mol), PDDG/MNDO (6.6 kcal/mol), PM3 (8.1 kcal/mol), MNDO' (8.5 kcal/mol), AM1 (11.1 kcal/mol), and MNDO (14.0 kcal/mol). For normal-valent halogen-containing molecules, the PDDG methods also provide improved heats of formation over MNDO/d. Hypervalent compounds were not included in the training set and improvements over the standard NDDO methods with sp basis sets were not obtained. For small haloalkanes, the PDDG methods yield more accurate heats of formation than are obtained from density functional theory (DFT) with the B3LYP and B3PW91 functionals using large basis sets. PDDG/PM3 and PM3' also give improved binding energies over the standard NDDO methods for complexes involving halide anions, and they are competitive with B3LYP/6-311++G(d,p) results including thermal corrections. Among the semiempirical methods studied, PDDG/PM3 also generates the best agreement with high-level ab initio G2 and CCSD(T) intrinsic activation energies for S(N)2 reactions involving methyl halides and halide anions. Finally, the MAEs in ionization potentials, dipole moments, and molecular geometries show that the parameter sets for the PDDG and reoptimized NDDO methods reduce the MAEs in heats of formation without compromising the other important QM observables.

A QM-MM interface between CHARMM and TURBOMOLE: implementation and application to systems in bulk phase and biologically active systems

The implementation of a hybrid QM-MM approach combining ab initio and density functional methods of TURBOMOLE with the molecular mechanics program package CHARMM is described. An interface has been created to allow data exchange between the two applications. With this method the efficient multiprocessor capabilities of TURBOMOLE can be utilized with CHARMM running as a single processor application. Therefore, features of nonparallel running code in CHARMM like the TRAVEL module for locating saddle points or VIBRAN for the calculation of second derivatives can be exploited by running the CPU intensive QM calculations in parallel. To test the methodology, several small systems are studied with both Hartree-Fock and density functional methods and varying QM-MM boundaries. Also, the computationally efficient RI-J method has been examined for use in QM-MM applications. A B(12) cofactor containing cobalt has been studied, to examine systems with a large QM region and transition metals. All tested methods perform satisfactory in comparison with pure quantum calculations. Additionally, algorithms for the characterization of saddle points have been tested for their potential use in QM-MM problems. The TRAVEL module of CHARMM has been applied to the Menshutkin reaction in the condensed phase, and a saddle point was located. This saddle point was verified by calculation of a steepest descent path connecting educt, transition state, and product, and by calculation of vibrational modes.

Some critical issues in the application of quantum mechanics/molecular mechanics methods to the study of transition metal complexes

The application of quantum mechanics/molecular mechanics (QM/MM) methods in transition metal chemistry is growing steadily. It becomes therefore appropriate to assess the importance of a number of technical issues associated to their implementation. This work presents the discussion of several of these issues, including the eventual need for conformational searches, the choice of the MM force field and the possibility of its tuning. The examples presented here prove that a proper handling of these technical aspects can lead to an improvement in the efficiency and quality of QM/MM calculations.

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.

New developments in applying quantum mechanics to proteins

Algorithmic improvements of quantum mechanical methodologies have increased our ability to study the electronic structure of fragments of a biomolecule (e.g. an enzyme active site) or entire biomolecules. Three main strategies have emerged as ways in which quantum mechanics can be applied to biomolecules. The supermolecule approach continues to be utilized, but it is slowly being replaced by the so-called coupled quantum mechanical/molecular mechanical methodologies. An exciting new direction is the continued development and application of linear-scaling quantum mechanical approaches to biomolecular systems.