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.     

Flexible-boundary QM/MM calculations: II. Partial charge transfer across the QM/MM boundary that passes through a covalent bond

Recently, based on the principle of electronic chemical potential equalization and the principle of charge conservation, we proposed a flexible-boundary scheme that allows both partial charge transfer and self-consistent polarization between the quantum mechanical (QM) and molecular mechanical (MM) subsystems in QM/MM calculations; the scheme was applied to study the atomic charges in selected ion–solvent complexes. In the present contribution, we further extend the flexible-boundary treatment to handle the QM/MM boundary passing through covalent bonds. We find that the flexible-boundary redistributed charge and dipole schemes yield reasonable agreement with full-QM calculations for a number of molecular ions and amino acids with charged side chains. Using the full-QM results as reference, the mean unsigned deviations are computed to be 0.06 e for atomic partial charges of the QM atoms, 0.11 e for the amounts of charge transfer between the QM and MM subsystems, and 0.016 Å for the lengths of the covalent bonds that directly connect the QM and MM subsystems. The results indicate the importance of accounting for partial charge transfer across the QM/MM boundary when the QM subsystems are charged.     

YinYang atom: a simple combined ab initio quantum mechanical molecular mechanical model

A simple interface is proposed for combined quantum mechanical (QM) molecular mechanical (MM) calculations for the systems where the QM and MM regions are connected through covalent bonds. Within this model, the atom that connects the two regions, called YinYang atom here, serves as an ordinary MM atom to other MM atoms and as a hydrogen-like atom to other QM atoms. Only one new empirical parameter is introduced to adjust the length of the connecting bond and is calibrated with the molecule propanol. This model is tested with the computation of equilibrium geometries and protonation energies for dozens of molecules. Special attention is paid on the influence of MM point charges on optimized geometry and protonation energy, and it is found that it is important to maintain local charge-neutrality in the MM region in order for the accurate calculation of the protonation and deprotonation energies. Overall the simple YinYang atom model yields comparable results to some other QM/MM models.     

A combined quantum mechanical and molecular mechanical method using modified generalized hybrid orbitals: implementation for electronic excited states

The generalized hybrid orbital (GHO) method is implemented at the second-order approximate coupled cluster singles and doubles (CC2) level for quantum mechanical (QM)/molecular mechanical (MM) electronic excited state calculations. The linear response function of CC2 in the GHO scheme is derived and implemented. The new implementation is applied to the first singlet excited states of three aromatic amino acids, phenylalanine, tyrosine, and tryptophan, and also bacteriorhodopsin for assessment. The results obtained for aromatic amino acids agreed well with the full QM CC2 calculations, while the calculated excitation energies of bacteriorhodopsin and its chromophore, all-trans retinal, reproduced the environmental shift of the experimental data. For the bacteriorhodopsin case, the environmental shift of GHO also showed good agreements with the experimental data. The contribution of the quantum effect of certain moieties in the excited states is elucidated by changing the partitioning of QM and MM regions.     

A critical evaluation of different QM/MM frontier treatments with SCC-DFTB as the QM method

The performance of different link atom based frontier treatments in QM/MM simulations was evaluated critically with SCC-DFTB as the QM method. In addition to the analysis of gas-phase molecules as in previous studies, an important element of the present work is that chemical reactions in realistic enzyme systems were also examined. The schemes tested include all options available in the program CHARMM for SCC-DFTB/MM simulation, which treat electrostatic interactions due to the MM atoms close to the QM/MM boundary in different ways. In addition, a new approach, the divided frontier charge (DIV), has been implemented in which the partial charge associated with the frontier MM atom ("link host") is evenly distributed to the other MM atoms in the same group. The performance of these schemes was evaluated based on properties including proton affinities, deprotonation energies, dipole moments, and energetics of proton transfer reactions. Similar to previous work, it was found that calculated proton affinities and deprotonation energies of alcohols, carbonic acids, amino acids, and model DNA bases are very sensitive to the link atom scheme; the commonly used single link atom approach often gives error on the order of 15 to 20 kcal/mol. Other schemes give better and, on average, mutually comparable results. For proton transfer reactions, encouragingly, both activation barriers and reaction energies are fairly insensitive (within a typical range of 2-4 kcal/mol) to the link atom scheme due to error cancellation, and this was observed for both gas-phase and enzyme systems. Therefore, the effect of using different link atom schemes in QM/MM simulations is rather small for chemical reactions that conserve the total charge. Although the current study used an approximate DFT method as the QM level, the observed trends are expected to be applicable to QM/MM methods with use of other QM approaches. This observation does not mean to encourage QM/MM simulations without careful benchmark in the study of specific systems, rather it emphasizes that other technical details, such as the treatment of long-range electrostatics, tend to play a more important role and need to be handled carefully.     

QM/MM calculation of protein magnetic shielding tensors with generalized hybrid-orbital method: a GIAO approach

A method to compute magnetic shielding tensors with generalized hybrid-orbital (GHO) QM/MM scheme is developed at the levels of Hartree-Fock and second-order M?ller-Plesset perturbation theory using gauge-including atomic orbitals. A feature of the GHO method is utilized to ensure gauge-origin independency of GHO shielding tensors in a simple way. The benchmark calculations indicate that the GHO method reproduced full-QM shielding constants nearly quantitatively for atoms not directly coupled to the GHO linking atoms. As an application to a realistic protein, carbon chemical shifts are calculated for the retinal chromophore in visual rhodopsin.     

Toward a new approach for determination of solute's charge distribution to analyze interatomic electrostatic interactions in quantum mechanical/molecular mechanical simulations

We present an alternative approach to determine "density-dependent property"-derived charges for molecules in the condensed phase. In the case of a solution, it is essential to take into consideration the electron polarization of molecules in the active site of this system. The solute and solvent molecules in this site have to be described by a quantum mechanical technique and the others are allowed to be treated by a molecular mechanical method (QM/MM scheme). For calculations based on this scheme, using the forces and interaction energy as density-dependent property our charges from interaction energy and forces (CHIEF) approach can provide the atom-centered charges on the solute atoms. These charges reproduce well the electrostatic potentials around the solvent molecules and present properly the picture of the electron density of the QM subsystem in the solution system. Thus, the CHIEF charges can be considered as the atomic charges under the conditions of the QM/MM simulation, and then enable one to analyze electrostatic interactions between atoms in the QM and MM regions. This approach would give a view of the QM nuclei and electrons different from the conventional methods.     

Performance assessment of semiempirical molecular orbital methods in describing halogen bonding: quantum mechanical and quantum mechanical/molecular mechanical-molecular dynamics study

The performance of semiempirical molecular-orbital methods--MNDO, MNDO-d, AM1, RM1, PM3 and PM6--in describing halogen bonding was evaluated, and the results were compared with molecular mechanical (MM) and quantum mechanical (QM) data. Three types of performance were assessed: (1) geometrical optimizations and binding energy calculations for 27 halogen-containing molecules complexed with various Lewis bases (Two of the tested methods, AM1 and RM1, gave results that agree with the QM data.); (2) charge distribution calculations for halobenzene molecules, determined by calculating the solvation free energies of the molecules relative to benzene in explicit and implicit generalized Born (GB) solvents (None of the methods gave results that agree with the experimental data.); and (3) appropriateness of the semiempirical methods in the hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme, investigated by studying the molecular inhibition of CK2 protein by eight halobenzimidazole and -benzotriazole derivatives using hybrid QM/MM molecular-dynamics (MD) simulations with the inhibitor described at the QM level by the AM1 method and the rest of the system described at the MM level. The pure MM approach with inclusion of an extra point of positive charge on the halogen atom approach gave better results than the hybrid QM/MM approach involving the AM1 method. Also, in comparison with the pure MM-GBSA (generalized Born surface area) binding energies and experimental data, the calculated QM/MM-GBSA binding energies of the inhibitors were improved by replacing the G(GB,QM/MM) solvation term with the corresponding G(GB,MM) term.     

Design-atom approach for the quantum mechanical/molecular mechanical covalent boundary: a design-carbon atom with five valence electrons

A critical issue underlying the accuracy and applicability of the combined quantum mechanical/molecular mechanical (QM/MM) methods is how to describe the QM/MM boundary across covalent bonds. Inspired by the ab initio pseudopotential theory, here we introduce a novel design atom approach for a more fundamental and transparent treatment of this QM/MM covalent boundary problem. The main idea is to replace the boundary atom of the active part with a design atom, which has a different number of valence electrons but very similar atomic properties. By modifying the Troullier-Martins scheme, which has been widely employed to construct norm-conserving pseudopotentials for density functional calculations, we have successfully developed a design-carbon atom with five valence electrons. Tests on a series of molecules yield very good structural and energetic results and indicate its transferability in describing a variety of chemical bonds, including double and triple bonds.     

Quantum mechanics/molecular mechanics electrostatic embedding with continuous and discrete functions

A quantum mechanics/molecular mechanics (QM/MM) implementation that uses the Gaussian electrostatic model (GEM) as the MM force field is presented. GEM relies on the reproduction of electronic density by using auxiliary basis sets to calculate each component of the intermolecular interaction. This hybrid method has been used, along with a conventional QM/MM (point charges) method, to determine the polarization on the QM subsystem by the MM environment in QM/MM calculations on 10 individual H(2)O dimers and a Mg(2+)-H(2)O dimer. We observe that GEM gives the correct polarization response in cases when the MM fragment has a small charge, while the point charges produce significant over-polarization of the QM subsystem and in several cases present an opposite sign for the polarization contribution. In the case when a large charge is located in the MM subsystem, for example, the Mg(2+) ion, the opposite is observed at small distances. However, this is overcome by the use of a damped Hermite charge, which provides the correct polarization response.     

Generalized hybrid orbital for the treatment of boundary atoms in combined quantum mechanical and molecular mechanical calculations using the semiempirical parameterized model 3 method

The application of combined quantum mechanical (QM) and molecular mechanical methods to large molecular systems requires an adequate treatment of the boundary between the two approaches. In this article, we extend the generalized hybrid orbital (GHO) method to the semiempirical parameterized model 3 (PM3) Hamiltonian combined with the CHARMM force field. The GHO method makes use of four hybrid orbitals, one of which is included in the QM region in self-consistent field optimization and three are treated as auxiliary orbitals that do not participate in the QM optimization, but they provide an effective electric field for interactions. An important feature of the GHO method is that the semiempirical parameters for the boundary atom are transferable, and these parameters have been developed for a carbon boundary atom consistent with the PM3 model. The combined GHO-PM3/CHARMM model has been tested on molecular geometry and proton affinity for a series of organic compounds.Acknowledgement We thank the National Institutes of Health for support of this research.Contribution to the Jacopo Tomasi Honorary Issue     

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.     

Importance of accurate charges in molecular docking: quantum mechanical/molecular mechanical (QM/MM) approach

The extent to which accuracy of electric charges plays a role in protein-ligand docking is investigated through development of a docking algorithm, which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations. In this algorithm, fixed charges of ligands obtained from force field parameterization are replaced by QM/MM calculations in the protein environment, treating only the ligands as the quantum region. The algorithm is tested on a set of 40 cocrystallized structures taken from the Protein Data Bank (PDB) and provides strong evidence that use of nonfixed charges is important. An algorithm, dubbed "Survival of the Fittest" (SOF) algorithm, is implemented to incorporate QM/MM charge calculations without any prior knowledge of native structures of the complexes. Using an iterative protocol, this algorithm is able in many cases to converge to a nativelike structure in systems where redocking of the ligand using a standard fixed charge force field exhibits nontrivial errors. The results demonstrate that polarization effects can play a significant role in determining the structures of protein-ligand complexes, and provide a promising start towards the development of more accurate docking methods for lead optimization applications.     

A regularized and renormalized electrostatic coupling Hamiltonian for hybrid quantum-mechanical-molecular-mechanical calculations

We describe a regularized and renormalized electrostatic coupling Hamiltonian for hybrid quantum-mechanical (QM)-molecular-mechanical (MM) calculations. To remedy the nonphysical QM/MM Coulomb interaction at short distances arising from a point electrostatic potential (ESP) charge of the MM atom and also to accommodate the effect of polarized MM atom in the coupling Hamiltonian, we propose a partial-wave expansion of the ESP charge and describe the effect of a s-wave expansion, extended over the covalent radius r(c), of the MM atom. The resulting potential describes that, at short distances, large scale cancellation of Coulomb interaction arises intrinsically from the localized expansion of the MM point charge and the potential self-consistently reduces to 1r(c) at zero distance providing a renormalization to the Coulomb energy near interatomic separations. Employing this renormalized Hamiltonian, we developed an interface between the Car-Parrinello molecular-dynamics program and the classical molecular-dynamics simulation program Groningen machine for chemical simulations. With this hybrid code we performed QM/MM calculations on water dimer, imidazole carbon monoxide (CO) complex, and imidazole-heme-CO complex with CO interacting with another imidazole. The QM/MM results are in excellent agreement with experimental data for the geometry of these complexes and other computational data found in literature.     

The generalized hybrid orbital method for combined quantum mechanical/molecular mechanical calculations: formulation and tests of the analytical derivatives

 Hybrid quantum mechanical (QM) and molecular mechanical (MM) potentials are becoming increasingly important for studying condensed-phase systems but one of the outstanding problems in the field has been how to treat covalent bonds between atoms of the QM and MM regions. Recently, we presented a generalized hybrid orbital (GHO) method that was designed to tackle this problem for hybrid potentials using semiempirical QM methods [Gao et al. (1998) J Phys Chem A 102: 4714–4721]. We tested the method on some small molecules and showed that it performed well when compared to the purely QM or MM potentials. In this article, we describe the formalism for the determination of the GHO energy derivatives and then present the results of more tests aimed at validating the model. These tests, involving the calculation of the proton affinities of some model compounds and a molecular dynamics simulation of a protein, indicate that the GHO method will prove useful for the application of hybrid potentials to solution-phase macromolecular systems. Received: 4 October 1999 / Accepted: 18 December 1999 / Published online: 5 June 2000     

The implementation of a fast and accurate QM/MM potential method in Amber

Version 9 of the Amber simulation programs includes a new semi-empirical hybrid QM/MM functionality. This includes support for implicit solvent (generalized Born) and for periodic explicit solvent simulations using a newly developed QM/MM implementation of the particle mesh Ewald (PME) method. The code provides sufficiently accurate gradients to run constant energy QM/MM MD simulations for many nanoseconds. The link atom approach used for treating the QM/MM boundary shows improved performance, and the user interface has been rewritten to bring the format into line with classical MD simulations. Support is provided for the PM3, PDDG/PM3, PM3CARB1, AM1, MNDO, and PDDG/MNDO semi-empirical Hamiltonians as well as the self-consistent charge density functional tight binding (SCC-DFTB) method. Performance has been improved to the point where using QM/MM, for a QM system of 71 atoms within an explicitly solvated protein using periodic boundaries and PME requires less than twice the cpu time of the corresponding classical simulation.