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).     

Dynamic structures of phosphodiesterase-5 active site by combined molecular dynamics simulations and hybrid quantum mechanical/molecular mechanical calculations

Various quantum mechanical/molecular mechanical (QM/MM) geometry optimizations starting from an x-ray crystal structure and from the snapshot structures of constrained molecular dynamics (MD) simulations have been performed to characterize two dynamically stable active site structures of phosphodiesterase-5 (PDE5) in solution. The only difference between the two PDE5 structures exists in the catalytic, second bridging ligand (BL2) which is HO- or H2O. It has been shown that, whereas BL2 (i.e. HO-) in the PDE5(BL2 = HO-) structure can really bridge the two positively charged metal ions (Zn2+ and Mg2+), BL2 (i.e. H2O) in the PDE5(BL2 = H2O) structure can only coordinate Mg2+. It has been demonstrated that the results of the QM/MM geometry optimizations are remarkably affected by the solvent water molecules, the dynamics of the protein environment, and the electronic embedding charges of the MM region in the QM part of the QMM/MM calculation. The PDE5(BL2 = H2O) geometries optimized by using the QM/MM method in different ways show strong couplings between these important factors. It is interesting to note that the PDE5(BL2 = HO-) and PDE5(BL2 = H2O) geometries determined by the QM/MM calculations neglecting these three factors are all consistent with the corresponding geometries determined by the QM/MM calculations that account for all of these three factors. These results suggest the overall effects of these three important factors on the optimized geometries can roughly cancel out. However, the QM/MM calculations that only account for some of these factors could lead to considerably different geometries. These results might be useful also in guiding future QM/MM geometry optimizations on other enzymes.     

QuanPol: A full spectrum and seamless QM/MM program

The quantum chemistry polarizable force field program (QuanPol) is implemented to perform combined quantum mechanical and molecular mechanical (QM/MM) calculations with induced dipole polarizable force fields and induced surface charge continuum solvation models. The QM methods include Hartree–Fock method, density functional theory method (DFT), generalized valence bond theory method, multiconfiguration self‐consistent field method, Møller–Plesset perturbation theory method, and time‐dependent DFT method. The induced dipoles of the MM atoms and the induced surface charges of the continuum solvation model are self‐consistently and variationally determined together with the QM wavefunction. The MM force field methods can be user specified, or a standard force field such as MMFF94, Chemistry at Harvard Molecular Mechanics (CHARMM), Assisted Model Building with Energy Refinement (AMBER), and Optimized Potentials for Liquid Simulations‐All Atom (OPLS‐AA). Analytic gradients for all of these methods are implemented so geometry optimization and molecular dynamics (MD) simulation can be performed. MD free energy perturbation and umbrella sampling methods are also implemented. © 2013 Wiley Periodicals, Inc.     

Geometry optimization with QM/MM,ONIOM, and other combined methods. I. Microiterations and constraints

Hybrid energy methods such as QM/MM and ONIOM, that combine different levels of theory into one calculation, have been very successful in describing large systems. Geometry optimization methods can take advantage of the partitioning of these calculations into a region treated at a quantum mechanical (QM) level of theory and the larger, remaining region treated by an inexpensive method such as molecular mechanics (MM). A series of microiterations can be employed to fully optimize the MM region for each optimization step in the QM region. Cartesian coordinates are used for the MM region and are chosen so that the internal coordinates of the QM region remain constant during the microiterations. The coordinates of the MM region are augmented to permit rigid body translation and rotation of the QM region. This is essential if any atoms in the MM region are constrained, but it also improves the efficiency of unconstrained optimizations. Because of the microiterations, special care is needed for the optimization step in the QM region so that the system remains in the same local valley during the course of the optimization. The optimization methodology with microiterations, constraints, and step-size control are illustrated by calculations on bacteriorhodopsin and other systems.     

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.     

Multistructural microiteration technique for geometry optimization and reaction path calculation in large systems

We propose a multistructural microiteration (MSM) method for geometry optimization and reaction path calculation in large systems. MSM is a simple extension of the geometrical microiteration technique. In conventional microiteration, the structure of the non‐reaction‐center (surrounding) part is optimized by fixing atoms in the reaction‐center part before displacements of the reaction‐center atoms. In this method, the surrounding part is described as the weighted sum of multiple surrounding structures that are independently optimized. Then, geometric displacements of the reaction‐center atoms are performed in the mean field generated by the weighted sum of the surrounding parts. MSM was combined with the QM/MM‐ONIOM method and applied to chemical reactions in aqueous solution or enzyme. In all three cases, MSM gave lower reaction energy profiles than the QM/MM‐ONIOM‐microiteration method over the entire reaction paths with comparable computational costs. © 2017 Wiley Periodicals, Inc.     

Cuby: An integrative framework for computational chemistry

Cuby is a computational chemistry framework written in the Ruby programming language. It provides unified access to a wide range of computational methods by interfacing external software and it implements various protocols that operate on their results. Using structured input files, elementary calculations can be combined into complex workflows. For users, Cuby provides a unified and userfriendly way to automate their work, seamlessly integrating calculations carried out in different computational chemistry programs. For example, the QM/MM module allows combining methods across the interfaced programs and the builtin molecular dynamics engine makes it possible to run a simulation on the resulting potential. For programmers, it provides high‐level, object‐oriented environment that allows rapid development and testing of new methods and computational protocols. The Cuby framework is available for download at http://cuby4.molecular.cz . © 2016 Wiley Periodicals, Inc.     

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.     

LICHEM: A QM/MM program for simulations with multipolar and polarizable force fields

We introduce an initial implementation of the LICHEM software package. LICHEM can interface with Gaussian, PSI4, NWChem, TINKER, and TINKER–HP to enable QM/MM calculations using multipolar/polarizable force fields. LICHEM extracts forces and energies from unmodified QM and MM software packages to perform geometry optimizations, single‐point energy calculations, or Monte Carlo simulations. When the QM and MM regions are connected by covalent bonds, the pseudo‐bond approach is employed to smoothly transition between the QM region and the polarizable force field. A series of water clusters and small peptides have been employed to test our initial implementation. The results obtained from these test systems show the capabilities of the new software and highlight the importance of including explicit polarization. © 2016 Wiley Periodicals, Inc.     

Fast generation of broken‐symmetry states in a large system including multiple iron–sulfur assemblies: Investigation of QM/MM energies,clusters charges,and spin populations

A density functional theory study is presented regarding the energetics and the Mulliken population analyses of a quantum mechanical/molecular mechanical (QM/MM) system including multiple iron–sulfur clusters in the QM region. The [FeFe]‐hydrogenase from Desulfovibrio desulfuricans was studied, and both the active site (an Fe₆S₆ assembly generally referred to as the H‐cluster) and an ancillary Fe₄S₄ site were treated at the BP86‐RI/TZVP level. The antiferromagnetic coupling that characterizes both sites was modeled using the broken‐symmetry (BS) approach. For such a QM system, 36 different BS couplings can be defined, depending on the localization of spin excess on the various spin centers. All the BS states were obtained by means of an effective and simple method for spin localization, that is here described and compared with more sophisticated approaches already available in literature. The variation of the QM/MM energy among the various geometry‐optimized protein models was found to be less than 25 kJ mol^–1. This energy variation almost doubles if no geometry optimization is performed. A detailed analysis of the additive nature of these variations in QM/MM energy is reported. The Mulliken charges show very small variations among the 36 BS states, whereas the Mulliken spin populations were found to be somewhat more variable. The relevance of such variations is discussed in light of the available Mössbauer and Electron Paramagnetic Resonance (EPR) spectroscopic data for the enzyme. Finally, the influence of the basis set on the spin populations, charges, and structural parameters of the models was investigated, by means of QM/MM computations on the same system at the BP86‐RI/SVP level. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

A QM/MM study on the aqueous solvation of the tetrahydroxouranylate [UO2(OH)4]2- complex ion

We report a QM augmented QM/MM study on the coordination of the tetrahydroxouranylate ion in aqueous solution. QM/MM geometry optimizations followed by full QM single-point calculations on the optimized structures show that a hexa-coordinated structure is more stable than the hepta-coordinated structure by 43 kJ/mol. Charge transfer of the tetrahydroxouranylate to the solvating water molecules is relatively modest, and can be modeled by including a solvation layer consisting of 12 explicit water molecules.     

Quantum capping potentials with point charges: a simple QM/MM approach for the calculation of large-molecule NMR shielding tensors

A simple quantum-mechanics/molecular-mechanics (QM/MM) approach for calculating NMR shielding tensors (sigma) is presented. The method involves capping the QM region with quantum capping potentials (QCPs) and representing the MM region with point charges. Test calculations on simple systems without MM charges show that calculated sigma values improve relative to the full QM results with increasing distance between the capped bond and chromophore. Calculations on the histidine amino acid and cytosine monophosphate (CMP) nucleic acid show that the use of QCPs with point charges result in mean errors in the isotropic component of sigma that are less than 1.6 ppm. The results also reveal that, contrary to previous work, the explicit effect of point charges on sigma through coupling with gauge factors, as in the gauge including atomic orbital approach, is minimal for the CMP molecule. The present QM/MM approach for calculating sigma is easy to apply and requires no code modification.     

Catalytic Mechanism of Cytochrome P450 2D6 for 4-Hydroxylation of Aripiprazole： A QM/MM Study

Drug metabolism is an important issue in drug discovery. Understanding how a drug is metabolized in the body will provide helpful information for lead optimization. Cytochrome P450 2D6 （CYP2D6） is a key enzyme for drug metabolism and responsible for the metabolism of about one third marketed drugs. Aripiprazole is an atypical an- tipsychotic and metabolized by CYP2D6 to its hydroxylated form. In this study, a series of computational methods were performed to understand how CYP2D6 accomplishes the 4-hydroxylation of aripiprazole. Molecular docking and molecular dynamics simulations were first performed to prepare the initial conformations for QM/MM calcula- tions. The results revealed two possible conformations for the drug-CYP2D6 complex. The ONIOM method for QM/MM calculations was then carried out to show detailed reaction pathways for the CYP2D6-catalyzed aripipra- zole hydroxylation reaction, which demonstrated that the dominant reactive channel was electrophilic and involved an initial attack on the n-system of the dichlorophenyl group of aripiprazole to produce cation δ-complex. Further- more, the product complex for each conformation was thermodynamically stable, which is in good agreement with previous reports.     

Contributions of pauli repulsions to the energetics and physical properties computed in QM/MM methods

Conventional combined quantum mechanical/molecular mechanical (QM/MM) methods lack explicit treatment of Pauli repulsions between the quantum‐mechanical and molecular‐mechanical subsystems. Instead, classical Lennard‐Jones (LJ) potentials between QM and MM nuclei are used to model electronic Pauli repulsion and long‐range London dispersion, despite the fact that the latter two are inherently of quantum nature. Use of the simple LJ potential in QM/MM methods can reproduce minimal geometries and energies of many molecular clusters reasonably well, as compared to full QM calculations. However, we show here that the LJ potential cannot correctly describe subtle details of the electron density of the QM subsystem because of the neglect of Pauli repulsions between the QM and MM subsystems. The inaccurate electron density subsequently affects the calculation of electronic and magnetic properties of the QM subsystem. To explicitly consider Pauli interactions with QM/MM methods, we propose a method to use empirical effective potentials on the MM atoms. The test case of the binding energy and magnetic properties of a water dimer shows promising results for the general application of effective potentials to mimic Pauli repulsions in QM/MM calculations. © 2013 Wiley Periodicals, Inc.     

QM:QM electronic embedding using Mulliken atomic charges: energies and analytic gradients in an ONIOM framework

An accurate first-principles treatment of chemical reactions for large systems remains a significant challenge facing electronic structure theory. Hybrid models, such as quantum mechanics:molecular mechanics (QM:MM) and quantum mechanics:quantum mechanics (QM:QM) schemes, provide a promising avenue for such studies. For many chemistries, including important reactions in materials science, molecular mechanics or semiempirical methods may not be appropriate, or parameters may not be available (e.g., surface chemistry of compound semiconductors such as indium phosphide or catalytic chemistry of transition metal oxides). In such cases, QM:QM schemes are of particular interest. In this work, a QM:QM electronic embedding model within the ONIOM (our own N-layer integrated molecular orbital molecular mechanics) extrapolation framework is presented. To define the embedding potential, we choose the real-system low-level Mulliken atomic charges. This results in a set of well-defined and unique embedding charges. However, the parametric dependence of the charges on molecular geometry complicates the energy gradient that is necessary for the efficient exploration of potential energy surfaces. We derive an efficient form for the forces where a single set of self-consistent field response equations is solved. Initial tests of the method and key algorithmic issues are discussed.     

Flexible effective fragment QM/MM method: validation through the challenging tests

A new version of the QM/MM method, which is based on the effective fragment potential (EFP) methodology [Gordon, M. et al., J Phys Chem A 2001, 105, 293] but allows flexible fragments, is verified through calculations of model molecular systems suggested by different authors as challenging tests for QM/MM approaches. For each example, the results of QM/MM calculations for a partitioned system are compared to the results of an all-electron ab initio quantum chemical study of the entire system. In each case we were able to achieve approximately similar or better accuracy of the QM/MM results compared to those described in original publications. In all calculations we kept the same set of parameters of our QM/MM scheme. A new test example is considered when calculating the potential of internal rotation in the histidine dipeptide around the C(alpha)bond;C(beta) side chain bond.