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.     

Structures of the Peptide–Water Complexes Studied by the Hybrid Quantum Mechanical—Molecular Mechanical (QM/MM) Technique

Applications of a new approach to the hybrid quantum mechanical and molecular mechanical (QM/MM) theory based on the effective fragment potential technique to calculations of the structures of the peptide—water complexes are described. Our approach assumes that the MM subsystem is viewed as a flexible composition of effective fragments, while fragment–fragment interactions are replaced by MM force fields. In this work, the QM subsystem is composed of water molecules and the MM part refers to peptides. Different isomers of the hydrogen-bonded complex of the dipeptide N-acetyl-L-alanine N-methylamide (AAMA) with four water molecules are considered, and the results of QM/MM calculations are compared to experimental data and to the results of the density functional theory (DFT) treatment. The properties of water chains inside polypeptide tubes, modeling proton wires inside ionic channels, are described.     

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.     

An extensible interface for QM/MM molecular dynamics simulations with AMBER

We present an extensible interface between the AMBER molecular dynamics (MD) software package and electronic structure software packages for quantum mechanical (QM) and mixed QM and classical molecular mechanical (MM) MD simulations within both mechanical and electronic embedding schemes. With this interface, ab initio wave function theory and density functional theory methods, as available in the supported electronic structure software packages, become available for QM/MM MD simulations with AMBER. The interface has been written in a modular fashion that allows straight forward extensions to support additional QM software packages and can easily be ported to other MD software. Data exchange between the MD and QM software is implemented by means of files and system calls or the message passing interface standard. Based on extensive tests, default settings for the supported QM packages are provided such that energy is conserved for typical QM/MM MD simulations in the microcanonical ensemble. Results for the free energy of binding of calcium ions to aspartate in aqueous solution comparing semiempirical and density functional Hamiltonians are shown to demonstrate features of this interface. © 2013 Wiley Periodicals, Inc.     

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.     

Multipole moments for embedding potentials: Exploring different atomic allocation algorithms

Polarizable quantum mechanical (QM)/molecular mechanics (MM)‐embedding methods are currently among the most promising methods for computationally feasible, yet reliable, production calculations of localized excitations and molecular response properties of large molecular complexes, such as proteins and RNA/DNA, and of molecules in solution. Our aim is to develop a computational methodology for distributed multipole moments and their associated multipole polarizabilities which is accurate, computationally efficient, and with smooth convergence with respect to multipole order. As the first step toward this goal, we herein investigate different ways of obtaining distributed atom‐centered multipole moments that are used in the construction of the electrostatic part of the embedding potential. Our objective is methods that not only are accurate and computationally efficient, but which can be consistently extended with site polarizabilities including internal charge transfer terms. We present a new way of dealing with well‐known problems in relation to the use of basis sets with diffuse functions in conventional atomic allocation algorithms, avoiding numerical integration schemes. Using this approach, we show that the classical embedding potential can be systematically improved, also when using basis sets with diffuse functions, and that very accurate embedding potentials suitable for QM/MM embedding calculations can be acquired. © 2016 Wiley Periodicals, Inc.     

Interfacing the GROMOS (bio)molecular simulation software to quantum‐chemical program packages

The newly implemented quantum‐chemical/molecular‐mechanical (QM/MM) functionality of the Groningen molecular simulation (GROMOS) software for (bio)molecular simulation is described. The implementation scheme is based on direct coupling of the GROMOS C++ software to executables of the quantum‐chemical program packages MNDO and TURBOMOLE, allowing for an independent further development of these packages. The new functions are validated for different test systems using program and model testing techniques. The effect of truncating the QM/MM electrostatic interactions at various QM/MM cutoff radii is discussed and the application of semiempirical versus density‐functional Hamiltonians for a solute molecule in aqueous solution is compared. © 2012 Wiley Periodicals, Inc.     

Development of massive multilevel molecular dynamics simulation program,platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions

A massively parallel program for quantum mechanical‐molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc‐pVDZ and B3LYP/cc‐pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6‐31G** calculations. We also performed excited QM/MM‐MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH‐insensitive and photo‐stable ultramarine fluorescent protein. Platypus accelerated on‐the‐fly excited‐state QM/MM‐MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50‐ps (200,000‐step) on‐the‐fly excited‐state QM/MM‐MD simulations for the SIRIUS in water. © 2016 The Authors. Journal of Computational Chemistry Published by 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.     

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.     

QM/MM study of aqueous solvation of the uranyl fluoride [UO2F4(2-)] complex

The aqueous solvation of the uranylfluoride complex [UO(2)F(4) (2-)] was studied using full quantum mechanical (QM) and hybrid QM/molecular mechanics (MM) methods. Inclusion of a complete first solvation shell was found necessary to reproduce the experimentally observed heptacoordination of uranium. An efficient and accurate computational model is proposed that consists of structure optimization of the coordinated uranium complex as QM region, followed by single-point full QM calculations to compute relative energies. This method is proven feasible for studies of large solvated actinide complexes.     

CPMD/GULP QM/MM interface for modeling periodic solids: Implementation and its application in the study of Y‐zeolite supported Rhn clusters

We report here the development of hybrid quantum mechanics/molecular mechanics (QM/MM) interface between the plane‐wave density functional theory based CPMD code and the empirical force‐field based GULP code for modeling periodic solids and surfaces. The hybrid QM/MM interface is based on the electrostatic coupling between QM and MM regions. The interface is designed for carrying out full relaxation of all the QM and MM atoms during geometry optimizations and molecular dynamics simulations, including the boundary atoms. Both Born–Oppenheimer and Car–Parrinello molecular dynamics schemes are enabled for the QM part during the QM/MM calculations. This interface has the advantage of parallelization of both the programs such that the QM and MM force evaluations can be carried out in parallel to model large systems. The interface program is first validated for total energy conservation and parallel scaling performance is benchmarked. Oxygen vacancy in α‐cristobalite is then studied in detail and the results are compared with a fully QM calculation and experimental data. Subsequently, we use our implementation to investigate the structure of rhodium cluster (Rh_n; n = 2 to 6) formed from Rh(C₂H₄)₂ complex adsorbed within a cavity of Y‐zeolite in a reducible atmosphere of H₂ gas. © 2016 Wiley Periodicals, Inc.     

Incorporation of deMon2k as a new parallel quantum mechanical code for the PUPIL system

The PUPIL system is a combination of software and protocols for the systematic linkage and interoperation of molecular dynamics and quantum mechanics codes to perform QM/MD (sometimes called QM/MM) calculations. The Gaussian03 and Amber packages were added to the PUPIL suite recently. However, efficient parallel QM codes are critical because calculation of the QM forces is the overwhelming majority of the computational load. Here we report details of incorporation of the deMon2k density functional suite as a new parallel QM code. An additional motivation is to add a highly optimized, purely DFT code. We illustrate with a demonstration study of the influence of perchlorate as a dopant ion of the poly(3,4‐ethylenedioxythiophene) conducting polymer in explicit acetonitrile solvent using Amber and deMon2k. We discuss unanticipated requirements for use of a scheme for semi‐empirical correction of Kohn‐Sham eigenvalues to give physically meaningful one‐electron gap energies. We provide comparison of both geometric parameters and electronic properties for nondoped and doped systems. We also present results comparing deMon2k and Gaussian03 calculation of forces for a short sequence of steps. We discuss briefly some difficult problems of quantum zone SCF convergence for the anionically doped system. The difficulties seem to be caused by well‐know deficiencies in simple approximate exchange‐correlation functionals. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

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.     

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.     

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.     

QM/MM calculation of solvent effects on absorption spectra of guanine

Electronic spectra of guanine in the gas phase and in water were studied by quantum mechanical/molecular mechanical (QM/MM) methods. Geometries for the excited‐state calculations were extracted from ground‐state molecular dynamics (MD) simulations using the self‐consistent‐charge density functional tight binding (SCC‐DFTB) method for the QM region and the TIP3P force field for the water environment. Theoretical absorption spectra were generated from excitation energies and oscillator strengths calculated for 50 to 500 MD snapshots of guanine in the gas phase (QM) and in solution (QM/MM). The excited‐state calculations used time‐dependent density functional theory (TDDFT) and the DFT‐based multireference configuration interaction (DFT/MRCI) method of Grimme and Waletzke, in combination with two basis sets. Our investigation covered keto‐N7H and keto‐N9H guanine, with particular focus on solvent effects in the low‐energy spectrum of the keto‐N9H tautomer. When compared with the vertical excitation energies of gas‐phase guanine at the optimized DFT (B3LYP/TZVP) geometry, the maxima in the computed solution spectra are shifted by several tenths of an eV. Three effects contribute: the use of SCC‐DFTB‐based rather than B3LYP‐based geometries in the MD snapshots (red shift of ca. 0.1 eV), explicit inclusion of nuclear motion through the MD snapshots (red shift of ca. 0.1 eV), and intrinsic solvent effects (differences in the absorption maxima in the computed gas‐phase and solution spectra, typically ca. 0.1–0.3 eV). A detailed analysis of the results indicates that the intrinsic solvent effects arise both from solvent‐induced structural changes and from electrostatic solute–solvent interactions, the latter being dominant. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010