Computing the (1)H NMR spectrum of a bulk ionic liquid from snapshots of car-parrinello molecular dynamics simulations.

We have investigated the performance of several computational protocols in predicting the NMR spectrum of a molecular ion in a complex liquid phase such as an ionic liquid. To do this, we computed the proton NMR chemical shifts of the 1-ethyl-3-methylimidazolium cation [emim](+) in [emim][Cl]. Environmental effects on the imidazolium ring proton chemical shifts are quite significant and must be taken into account explicitly. Calculations performed on the isolated imidazolium cation as well as on the [emim][Cl] ion pair grossly fail to reproduce the correct spacing between proton signals. In contrast, calculations performed on clusters extracted from the trajectory of a Car-Parrinello molecular dynamics simulation yield very good results.     

A new QM/MM method oriented to the study of ionic liquids

The interest on room temperature ionic liquids has grown in the last decades because of their use as all‐purpose solvent and their low environmental impact. In the present work, a new theoretical procedure is developed to study pure ionic liquids within the framework of the quantum mechanics/molecular mechanics method. Each type of ion (cation or anion) is considered as an independent entity quantum mechanically described that follows a differentiated path in the liquid. The method permits, through an iterative procedure, the full coupling between the polarized charge distribution of the ions and the liquid structure around them. The procedure has been tested with 1‐ethyl‐3‐methylimidazolium tetrafluoroborate. It was found that, similar to non‐polar liquids and as a consequence of the low value of the reaction field, the cation and anion charge distributions are hardly polarized by the rest of molecules in the liquid. Their structure is characterized by an alternance between anion and cation shells as evidenced by the coincidence of the first maximum of the anion–anion and cation–cation radial distribution functions with the first minimum of the anion‐cation. Some degree of stacking between the cations is also found. © 2015 Wiley Periodicals, Inc.     

Prediction and rationalization of protein pKa values using QM and QM/MM methods

We describe the development and application of a computational method for the prediction and rationalization of pKa values of ionizable residues in proteins, based on ab initio quantum mechanics (QM) and the effective fragment potential (EFPs) method (a hybrid QM/MM method). The theoretical developments include (1) a covalent boundary method based on frozen localized orbitals, (2) divide-and-conquer methods for the ab initio computation of protein EFPs consisting of multipoles up to octupoles and dipole polarizability tensors, (3) a method for computing vibrational free energies for a localized molecular region, and (4) solutions of the polarized continuum model of bulk solvation equations for protein-sized systems. The QM-based pKa prediction method is one of the most accurate methods currently available and can be used in cases where other pKa prediction methods fail. Preliminary analysis of the computed results indicate that many pKa values (1) are primarily determined by hydrogen bonds rather than long-range charge-charge interactions and (2) are relatively insensitive to large-scale dynamical fluctuations of the protein structure.     

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.     

Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model

The most general way to improve the accuracy of binding‐affinity calculations for protein–ligand systems is to use quantum‐mechanical (QM) methods together with rigorous alchemical‐perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin‐3 at a reasonably high QM level (dispersion‐corrected density functional theory with a triple‐zeta basis set) and with a sufficiently large QM system to include all short‐range interactions with the ligand (744–748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non‐Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ～3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations. © 2015 Wiley Periodicals, Inc.     

Hybrid QM/MM studies on energetics of malonaldehyde in condensed phase

Energetics of eight enol isomers of a malonaldehyde (MA) at the ground state in CCl₄ solvent environment have been investigated using a hybrid quantum mechanical (QM)/molecular mechanics (MM) method. It is found that relative energies of the isomers slightly change due to interactions between MA and surrounding atoms. In an isolated environment all eight isomers have stable planar structure. On the other hand, most of the isomers have nonplanar structure in CCl₄ whose interactions with the solute molecules are, however, weak. Mainly, structural changes are found in the H? O? C?C dihedral angle ?_{H? O? C?C}, i.e., a hydrogen atom that is connected directly with an oxygen atom, is located at a nonplanar position, and other atoms remain almost planar. Vertical excitation energies of low‐lying excited states at the resultant optimized structure of each isomer are evaluated. The vertical excitation energies in CCl₄ are almost the same as those in the isolated environment, but some changes were found in the triplet excitation states. © 2004 Wiley Periodicals, Inc. Int J Quantum Chem, 2004     

Revealing the molecular mechanisms of proteolysis of SARS-CoV-2 Mpro by QM/MM computational methods

SARS-CoV-2 M^pro is one of the enzymes essential for the replication process of the virus responsible for the COVID-19 pandemic. This work is focused on exploring its proteolysis reaction by means of QM/MM methods. The resulting free energy landscape of the process provides valuable information on the species appearing along the reaction path and suggests that the mechanism of action of this enzyme, taking place in four steps, slightly differs from that of other cysteine proteases. Our predictions, which are in agreement with some recently published experimental data, can be used to guide the design of COVID-19 antiviral compounds with clinical potential.

The molecular mechanism of the proteolysis reaction catalyzed by SARS-CoV-2 M^pro, one of the enzymes essential for the replication process of the virus responsible for the COVID-19 pandemic, is described using computational QM/MM methods.     

51V NMR chemical shifts calculated from QM/MM models of vanadium chloroperoxidase

(51)V NMR chemical shifts calculated from QM/MM-optimized (QM=quantum mechanical; MM=molecular mechanical) models of vanadium-dependent chloroperoxidase (VCPO) are presented. An extensive number of protonation states for the vanadium cofactor (active site of the protein) and a number of probable positional isomers for each of the protonation states are considered. The size of the QM region is increased incrementally to observe the convergence behavior of the (51)V NMR chemical shifts. A total of 40 models are assessed by comparison to experimental solid-state (51)V NMR results recently reported in the literature. Isotropic chemical shifts are found to be a poor indicator of the protonation state; however, anisotropic chemical shifts and the nuclear quadrupole tensors appear to be sensitive to changes in the proton environment of the vanadium nuclei. This detailed investigation of the (51)V NMR chemical shifts computed from QM/MM models provides further evidence that the ground state is either a triply protonated (one axial water and one equatorial hydroxyl group) or a doubly protonated vanadate moiety in VCPO. Particular attention is given to the electrostatic and geometric effects of the protein environment. This is the first study to compute anisotropic NMR chemical shifts from QM/MM models of an active metalloprotein for direct comparison with solid-state MAS NMR data. This theoretical approach enhances the potential use of experimental solid-state NMR spectroscopy for the structural determination of metalloproteins.     

Chaperoned alchemical free energy simulations: a general method for QM, MM, and QM/MM potentials

A general method for alchemical free energy simulations using QM, MM, and QM/MM potential is developed by introducing "chaperones" to restrain the structures, particularly near the end points. A calculation of the free energy difference between two triazole tautomers in aqueous solution is used to illustrate the method.     

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.     

Reduction potentials and acidity constants of Mn superoxide dismutase calculated by QM/MM free-energy methods

We used two theoretical methods to estimate reduction potentials and acidity constants in Mn superoxide dismutase (MnSOD), namely combined quantum mechanical and molecular mechanics (QM/MM) thermodynamic cycle perturbation (QTCP) and the QM/MM-PBSA approach. In the latter, QM/MM energies are combined with continuum solvation energies calculated by solving the Poisson-Boltzmann equation (PB) or by the generalised Born approach (GB) and non-polar solvation energies calculated from the solvent-exposed surface area. We show that using the QTCP method, we can obtain accurate and precise estimates of the proton-coupled reduction potential for MnSOD, 0.30±0.01 V, which compares favourably with experimental estimates of 0.26-0.40 V. However, the calculated potentials depend strongly on the DFT functional used: The B3LYP functional gives 0.6 V more positive potentials than the PBE functional. The QM/MM-PBSA approach leads to somewhat too high reduction potentials for the coupled reaction and the results depend on the solvation model used. For reactions involving a change in the net charge of the metal site, the corresponding results differ by up to 1.3 V or 24 pK(a) units, rendering the QM/MM-PBSA method useless to determine absolute potentials. However, it may still be useful to estimate relative shifts, although the QTCP method is expected to be more accurate.     

QM/MM study of the taxadiene synthase mechanism

Combined quantum mechanics/molecular mechanics (QM/MM) calculations were used to investigate the reaction mechanism of taxadiene synthase (TXS). TXS catalyzes the cyclization of geranylgeranyl diphosphate (GGPP) to taxadiene (T) and four minor cyclic products. All these products originate from the deprotonation of carbocation intermediates. The reaction profiles for the conversion of GGPP to T as well as to minor products were calculated for different configurations of relevant TXS carbocation complexes. The QM region was treated at the M06-2X/TZVP level, while the CHARMM27 force field was used to describe the MM region. The QM/MM calculations suggest a reaction pathway for the conversion of GGPP to T, which slightly differs from previous proposals regarding the number of reaction steps and the conformation of the carbocations. The QM/MM results also indicate that the formation of minor products via water-assisted deprotonation of the carbocations is highly exothermic, by about −7 to −23 kcal/mol. Curiously, however, the computed barriers and reaction energies indicate that the formation of some of the minor products is more facile than the formation of T. Thus, the present QM/MM calculations provide detailed insights into possible reaction pathways and into the origin of the promiscuity of TXS, but they do not reproduce the product distribution observed experimentally. © 2019 Wiley Periodicals, Inc.     

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

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

Intermediates in dioxygen activation by methane monooxygenase: a QM/MM study

Protein effects in the activation of dioxygen by methane monooxygenase (MMO) were investigated by using combined QM/MM and broken-symmetry Density Functional Theory (DFT) methods. The effects of a novel empirical scheme recently developed by our group on the relative DFT energies of the various intermediates in the catalytic cycle are investigated. Inclusion of the protein leads to much better agreement between the experimental and computed geometric structures for the reduced form (MMOH(red)). Analysis of the electronic structure of MMOH(red) reveals that the two iron atoms have distinct environments. Different coordination geometries tested for the MMOH(peroxo) intermediate reveal that, in the protein environment, the mu-eta2,eta2 structure is more stable than the others. Our analysis also shows that the protein helps to drive reactants toward products along the reaction path. Furthermore, these results demonstrate the importance of including the protein environment in our models and the usefulness of the QM/MM approach for accurate modeling of enzymatic reactions. A discrepancy remains in our calculation of the Fe-Fe distance in our model of HQ as compared to EXAFS data obtained several years ago, for which we currently do not have an explanation.