QM/MM: what have we learned,where are we,and where do we go from here?

This paper briefly reviews the current status of the most popular methods for combined quantum mechanical/molecular mechanical (QM/MM) calculations, including their advantages and disadvantages. There is a special emphasis on very general link-atom methods and various ways to treat the charge near the boundary. Mechanical and electric embedding are contrasted. We consider methods applicable to gas-phase organic chemistry, liquid-phase organic and organometallic chemistry, biochemistry, and solid-state chemistry. Then we review some recent tests of QM/MM methods and summarize what we learn about QM/MM from these studies. We also discuss some available software. Finally, we present a few comments about future directions of research in this exciting area, where we focus on more intimate blends of QM with MM. Contribution to the Proceedings of the 10th Electronic Computational Chemistry Conference     

Combining crystallographic and quantum chemical data to understand DNA-protein π-interactions in nature

Noncovalent interactions are accepted to be prevalent across biochemical systems, including governing interactions between nucleic acids and proteins. The present review summarizes work done to characterize the abundance, structure and strength of DNA–protein π interactions by combining rigorous searches of experimental X-ray crystal structures of DNA–protein complexes and quantum chemical calculations. Focus is placed on interactions that occur between the π-containing amino acids (W, H, F, Y, R, E, and D) and the canonical DNA nucleobases (A, T, G, and C) or 2′-deoxyribose moiety. These studies highlight the considerable frequency of both DNA–protein π–π and sugar–π interactions in nature, which can involve any π-containing amino acid arranged in many unique binding orientations with respect to any DNA component. When combined with the significant strength predicted for the identified DNA–protein π contacts using density functional theory, these works underscore the potential impact of these interactions on critical biological functions. This conclusion is supported by a review of examples from the recent literature that have acknowledged the role of DNA–protein π interactions in binding, specificity, and catalysis.

     

Force field parameters for rotation around chi torsion axis in nucleic acids

To raise the accuracy of the force field for nucleic acids, several parameters were elaborated, focusing on the rotation around chi torsion axis. The reliability of molecular dynamics (MD) simulation was significantly increased by improving the torsion parameters at C8--N9--C1'--X (X = H1', C2', O4') in A, G and those at C6--N1--C1'--X in C, T, and U. In this work, we constructed small models representing the chemical structure of A, G, C, T, and U, and estimated energy profile for chi-axis rotation by executing numerous quantum mechanical (QM) calculations. The parameters were derived by discrete Fourier transformation of the calculated QM data. A comparison in energy profile between molecular mechanical (MM) calculation and QM one shows that our presently derived parameters well reproduce the energy surface of QM calculation for all the above torsion terms. Furthermore, our parameters show a good performance in MD simulations of some nucleic acids. Hence, the present refinement of parameters will enable us to perform more accurate simulations for various types of nucleic acids.     

OMx-D: semiempirical methods with orthogonalization and dispersion corrections. Implementation and biochemical application

The semiempirical methods of the OMx family (orthogonalization models OM1, OM2, and OM3) are known to describe biochemical systems more accurately than standard semiempirical approaches such as AM1. We investigate the benefits of augmenting these methods with an empirical dispersion term (OMx-D) taken from recent density functional work, without modifying the standard OMx parameters. Significant improvements are achieved for non-covalent interactions, with mean unsigned errors of 1.41 kcal/mol (OM2-D) and 1.31 kcal/mol (OM3-D) for the binding energy of the complexes in the JSCH-2005 data base. This supports the use of these augmented methods in quantum mechanical/molecular mechanical (QM/MM) studies of biomolecules, for example during system preparation and equilibration. As an illustrative application, we present QM and QM/MM calculations on the binding between antibody 34E4 and a hapten, where OM3-D performs better than the methods without dispersion terms (AM1, OM3).     

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.     

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.     

Density functional theory with dispersion corrections for supramolecular structures, aggregates, and complexes of (bio)organic molecules

Kohn-Sham density functional theory (KS-DFT) is nowadays the most widely used quantum chemical method for electronic structure calculations in chemistry and physics. Its further application in e.g. supramolecular chemistry or biochemistry has mainly been hampered by the inability of almost all current density functionals to describe the ubiquitous attractive long-range van der Waals (dispersion) interactions. We review here methods to overcome this defect, and describe in detail a very successful correction that is based on damped -C(6).R(-6) potentials (DFT-D). As examples we consider the non-covalent inter- and intra-molecular interactions in unsaturated organic molecules (so-called pi-pi stacking in benzenes and dyes), in biologically relevant systems (nucleic acid bases/pairs, proteins, and 'folding' models), between fluorinated molecules, between curved aromatics (corannulene and carbon nanotubes) and small molecules, and for the encapsulation of methane in water clusters. In selected cases we partition the interaction energies into the most relevant contributions from exchange-repulsion, electrostatics, and dispersion in order to provide qualitative insight into the binding character.     

Leading RNA tertiary interactions: structures, energies, and water insertion of A-minor and P-interactions. A quantum chemical view

Complex molecular shapes of ribosomal RNA molecules are stabilized by recurrent types of tertiary interactions involving highly specific and conserved non-Watson-Crick base pairs, triplets, and quartets. We analyzed the intrinsic structure and stability of the P-motif and the four basic types of A-minor interactions (types I, II, III, and 0), which represent the most prominent RNA tertiary interaction patterns refined in the course of evolution. In the studied interactions, the electron correlation component of the stabilization usually exceeds the Hartree-Fock (HF) term, leading to a strikingly different balance of forces as compared to standard base pairing stabilized primarily by the HF term. In other words, the A-minor and P-interactions are considerably more influenced by the dispersion energy as compared to canonical base pairs, which makes them particularly suitable to zip the folded RNA structures that are substantially hydrated even in their interior. Continuum solvent COSMO calculations confirm that the stability of the canonical GC base pair is affected (reduced) by the continuum solvent screening considerably more than the stability of the A-minor interaction. Among the studied systems, the strong A-minor II and weak A-minor III interactions require water molecules to stabilize the experimental geometry. Gas-phase optimization of the canonical A-minor II A/CG triplet without water results in a geometry that is clearly inconsistent with the RNA structure. The gas-phase structure of the P-interaction and the most stable A-minor I interaction nicely agrees with the geometries occurring in the ribosome. A-minor I can also adopt an alternative water-mediated substate rather often observed in X-ray and molecular dynamics studies. The A-minor I water bridge, however, does not appear to stabilize the tertiary contact, and its role is to provide structural flexibility to this binding pattern within the context of the RNA structure. Interestingly, the insertion of a polar water molecule in the A-minor I A/CG tertiary contact occurring in the A/C tertiary pair is stabilized primarily by the HF (electrostatic) interaction energy, while the dispersion-controlled A/G contact remains firmly bound. Thus, the intrinsic balance of forces as revealed by quantum mechanics (QM) calculations nicely correlates with many behavioral aspects of the studied interactions inside RNA. The comparison of interaction energies computed using quantum chemistry and an AMBER force field reveals that common molecular mechanics calculations perform rather well, except that the strength of the P-interaction is modestly overestimated. We also briefly discuss the non-negligible methodological differences when evaluating simple base-base nucleic acids base pairs and the complex RNA tertiary base pairing patterns using QM procedures.     

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.     

A QM/MM study of cisplatin-DNA oligonucleotides: from simple models to realistic systems

QM/MM calculations were employed to investigate the role of hydrogen bonding and pi stacking in several single- and double-stranded cisplatin-DNA structures. Computed geometrical parameters reproduce experimental structures of cisplatin and its complex with guanine-phosphate-guanine. Following QM/MM optimisation, single-point DFT calculations allowed estimation of intermolecular forces through atoms in molecules (AIM) analysis. Binding energies of platinated single-strand DNA qualitatively agree with myriad experimental and theoretical studies showing that complexes of guanine are stronger than those of adenine. The topology of all studied complexes confirms that platination strongly affects the stability of both single- and double-stranded DNAs: Pt-N-H...X (X = N or O) interactions are ubiquitous in these complexes and account for over 70 % of all H-bonding interactions. The pi stacking is greatly reduced by both mono- and bifunctional complexation: the former causes a loss of about 3-4 kcal mol(-1), whereas the latter leads to more drastic disruption. The effect of platination on Watson-Crick GC is similar to that found in previous studies: major redistribution of energy occurs, but the overall stability is barely affected. The BH&H/AMBER/AIM approach was also used to study platination of a double-stranded DNA octamer d(CCTG*G*TCC)d(GGACCAGG), for which an experimental structure is available. Comparison between theory and experiment is satisfactory, and also reproduces previous DFT-based studies of analogous structures. The effect of platination is similar to that seen in model systems, although the effect on GC pairing was more pronounced. These calculations also reveal weaker, secondary interactions of the form Pt...O and Pt...N, detected in several single- and double-stranded DNA.     

The effect of methylation on the hydrogen-bonding and stacking interaction of nucleic acid bases

Methylated nucleosides play an important role in DNA/RNA function, and may affect republication by interrupting the base-pairing and base-stacking. In order to investigate the effect of methylation on the interaction between nucleic acid bases, this work presents the hydrogen-bonding and stacking interactions between 5-methylcytosine and guanine (G), cytosine (C) and G, 1-methyladenine and thymine (T), as well as adenine and T. Geometry optimization and potential energy surface scan have been performed for the involved complexes by MP2 calculations. The interaction energies, which were corrected for the basis-set superposition error by the full Boys–Bernardi counterpoise correction scheme, were used to evaluate the interaction intensity of these nucleic acid bases. The atoms in molecules theory and natural bond orbital analysis have been performed to study the hydrogen bonds in these complexes. The result shows that the methyl substitute contributes the stability to these complexes because it enhances either the hydrogen bonding or the staking interaction between nucleic acid bases studied.     

Tautomerism of pyrimidines and purines in the gas phase and in low-temperature matrices,and some biological implications

Infrared spectra in the NH, OH, and C?O stretching regions are reviewed for various methylated and halogenated 2(4)-oxopyrimidines and uracils, as well as other nucleic acid derivatives, in the vapor phase and in low-temperature matrices. The 2-oxopyrimidines are predominantly in the enol form both in the vapor phase and in low-temperature matrices. By contrast, the 4-oxopyrimidines exhibit comparable populations of the keto and enol forms, with K_T ≈ 1–2, and a difference in chemical binding energy between the two forms in the gaseous phase of the order of 1–2 kcal/mol. The observed tautomeric equilibria in the gas phase point to the need for a drastic revision of interpretations of theoretical methods, and simultaneously provide the appropriate quantitative data necessary for testing the results of quantum mechanical calculations. In sharp contrast to other heterocyclic systems, several of the bases found in natural nucleic acids were found to exist predominantly in the keto or amino forms, as in the solution phase. In particular, uracils exist predominantly in the keto form. This has made possible, for this class of compounds, to evaluate the heats of sublimation and vaporization, and to relate these data to hydrogen bonding and stacking interactions in the condensed phases. Examples are presented of base analogs which do exhibit appreciable tautomerism in solution. Some biological implications of the foregoing are presented in relation to the types of heterocyclic bases found in natural nucleic acids, and to concepts of spontaneous and induced mutagenesis.     

Amino acid analogues bind to carbon nanotube via π-π interactions: comparison of molecular mechanical and quantum mechanical calculations

Understanding the interaction between carbon nanotubes (CNTs) and biomolecules is essential to the CNT-based nanotechnology and biotechnology. Some recent experiments have suggested that the π-π stacking interactions between protein's aromatic residues and CNTs might play a key role in their binding, which raises interest in large scale modeling of protein-CNT complexes and associated π-π interactions at atomic detail. However, there is concern on the accuracy of classical fixed-charge molecular force fields due to their classical treatments and lack of polarizability. Here, we study the binding of three aromatic residue analogues (mimicking phenylalanine, tyrosine, and tryptophan) and benzene to a single-walled CNT, and compare the molecular mechanical (MM) calculations using three popular fixed-charge force fields (OPLSAA, AMBER, and CHARMM), with quantum mechanical (QM) calculations using the density-functional tight-binding method with the inclusion of dispersion correction (DFTB-D). Two typical configurations commonly found in π-π interactions are used, one with the aromatic rings parallel to the CNT surface (flat), and the other perpendicular (edge). Our calculations reveal that compared to the QM results the MM approaches can appropriately reproduce the strength of π-π interactions for both configurations, and more importantly, the energy difference between them, indicating that the various contributions to π-π interactions have been implicitly included in the van der Waals parameters of the standard MM force fields. Meanwhile, these MM models are less accurate in predicting the exact structural binding patterns (matching surface), meaning there are still rooms to be improved. In addition, we have provided a comprehensive and reliable QM picture for the π-π interactions of aromatic molecules with CNTs in gas phase, which might be used as a benchmark for future force field developments.     

An application of the van der Waals density functional: Hydrogen bonding and stacking interactions between nucleobases

We apply the van der Waals density functional (vdW-DF) to study hydrogen bonding and stacking interactions between nucleobases. The excellent agreement of our results with high level quantum chemical calculations highlights the value of the vdW-DF for first-principles investigations of biologically important molecules. Our results suggest that, in the case of hydrogen-bonded nucleobase pairs, dispersion interactions reduce the cost of propeller twists while having a negligible effect on buckling. Furthermore, the efficient scaling of DFT methods allowed for the easy optimization of separation distance between nucleobase stacks, indicating enhancements in the interaction energy of up to 3 kcalmol over previous fixed distance calculations. We anticipate that these results are significant for extending the vdW-DF method to model larger vdW complexes and biological molecules.     

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.     

Computational and data driven molecular material design assisted by low scaling quantum mechanics calculations and machine learning

Electronic structure methods based on quantum mechanics (QM) are widely employed in the computational predictions of the molecular properties and optoelectronic properties of molecular materials. The computational costs of these QM methods, ranging from density functional theory (DFT) or time-dependent DFT (TDDFT) to wave-function theory (WFT), usually increase sharply with the system size, causing the curse of dimensionality and hindering the QM calculations for large sized systems such as long polymer oligomers and complex molecular aggregates. In such cases, in recent years low scaling QM methods and machine learning (ML) techniques have been adopted to reduce the computational costs and thus assist computational and data driven molecular material design. In this review, we illustrated low scaling ground-state and excited-state QM approaches and their applications to long oligomers, self-assembled supramolecular complexes, stimuli-responsive materials, mechanically interlocked molecules, and excited state processes in molecular aggregates. Variable electrostatic parameters were also introduced in the modified force fields with the polarization model. On the basis of QM computational or experimental datasets, several ML algorithms, including explainable models, deep learning, and on-line learning methods, have been employed to predict the molecular energies, forces, electronic structure properties, and optical or electrical properties of materials. It can be conceived that low scaling algorithms with periodic boundary conditions are expected to be further applicable to functional materials, perhaps in combination with machine learning to fast predict the lattice energy, crystal structures, and spectroscopic properties of periodic functional materials.

Low scaling quantum mechanics calculations and machine learning can be employed to efficiently predict the molecular energies, forces, and optical and electrical properties of molecular materials and their aggregates.