Quantum Mechanical Interaction Energies Research Articles

Accurate Reproduction of Quantum Mechanical Many-Body Interactions in Peptide Main-Chain Hydrogen-Bonding Oligomers by the Polarizable Gaussian Multipole Model.

A key advantage of polarizable force fields is their ability to model the atomic polarization effects that play key roles in the atomic many-body interactions. In this work, we assessed the accuracy of the recently developed polarizable Gaussian Multipole (pGM) models in reproducing quantum mechanical (QM) interaction energies, many-body interaction energies, as well as the nonadditive and additive contributions to the many-body interactions for peptide main-chain hydrogen-bonding conformers, using glycine dipeptide oligomers as the model systems. Two types of pGM models were considered, including that with (pGM-perm) and without (pGM-ind) permanent atomic dipoles. The performances of the pGM models were compared with several widely used force fields, including two polarizable (Amoeba13 and ff12pol) and three additive (ff19SB, ff15ipq, and ff03) force fields. Encouragingly, the pGM models outperform all other force fields in terms of reproducing QM interaction energies, many-body interaction energies, as well as the nonadditive and additive contributions to the many-body interactions, as measured by the root-mean-square errors (RMSEs) and mean absolute errors (MAEs). Furthermore, we tested the robustness of the pGM models against polarizability parameterization errors by employing alternative polarizabilities that are either scaled or obtained from other force fields. The results show that the pGM models with alternative polarizabilities exhibit improved accuracy in reproducing QM many-body interaction energies as well as the nonadditive and additive contributions compared with other polarizable force fields, suggesting that the pGM models are robust against the errors in polarizability parameterizations. This work shows that the pGM models are capable of accurately modeling polarization effects and have the potential to serve as templates for developing next-generation polarizable force fields for modeling various biological systems.

A Reduced Generalized Force Field for Biological Halogen Bonds.

The halogen bond (or X-bond) is a noncovalent interaction that is increasingly recognized as an important design tool for engineering protein-ligand interactions and controlling the structures of proteins and nucleic acids. In the past decade, there have been significant efforts to characterize the structure-energy relationships of this interaction in macromolecules. Progress in the computational modeling of X-bonds in biological molecules, however, has lagged behind these experimental studies, with most molecular mechanics/dynamics-based simulation methods not properly treating the properties of the X-bond. We had previously derived a force field for biological X-bonds (ffBXB) based on a set of potential energy functions that describe the anisotropic electrostatic and shape properties of halogens participating in X-bonds. Although fairly accurate for reproducing the energies within biomolecular systems, including X-bonds engineered into a DNA junction, the ffBXB with its seven variable parameters was considered to be too unwieldy for general applications. In the current study, we have generalized the ffBXB by reducing the number of variables to just one for each halogen type and show that this remaining electrostatic variable can be estimated for any new halogenated molecule through a standard restricted electrostatic potential calculation of atomic charges. In addition, we have generalized the ffBXB for both nucleic acids and proteins. As a proof of principle, we have parameterized this reduced and more general ffBXB against the AMBER force field. The resulting parameter set was shown to accurately recapitulate the quantum mechanical landscape and experimental interaction energies of X-bonds incorporated into DNA junction and T4 lysozyme model systems. Thus, this reduced and generalized ffBXB is more readily adaptable for incorporation into classical molecular mechanics/dynamics algorithms, including those commonly used to design inhibitors against therapeutic targets in medicinal chemistry and materials in biomolecular engineering.

Further Optimization and Validation of the Classical Drude Polarizable Protein Force Field.

The CHARMM Drude-2013 polarizable force field (FF) was developed to include the explicit treatment of induced electronic polarizability, resulting in a more accurate description of the electrostatic interactions in molecular dynamics (MD) simulations. While the Drude-2013 protein FF has shown success in improving the folding properties of α-helical peptides and to reproduce experimental observables in simulations up to 1 μs, some limitations were noted regarding the stability of β-sheet structures in simulations longer than 100 ns as well as larger deviations from crystal structures in simulations of a number of proteins compared to the additive CHARMM36 protein FF. The origin of the instability has been identified and appears to be primarily due to overestimated atomic polarizabilities and induced dipole-dipole interactions on the Cβ, Cγ, and Cδ side chain atoms. To resolve this and other issues, a number of aspects of the model were revisited, resulting in Drude-2019 protein FF. Backbone parameters were optimized targeting the conformational properties of the (Ala)5 peptide in solution along with gas phase properties of the alanine dipeptide. Dipeptides that contain N-acetylated and N'-methylamidated termini, excluding Gly, Pro, and Ala, were used as models to optimize the atomic polarizabilities and Thole screening factors on selected Cβ, Cγ, and Cδ carbons by targeting quantum mechanical (QM) dipole moments and molecular polarizabilities. In addition, to obtain better conformational properties, side chain χ1 and χ2 dihedral parameters were optimized targeting QM data for the respective side chain dipeptide conformations as well as Protein Data Bank survey data based on the χ1, χ2 sampling from Hamiltonian replica-exchange MD simulations of (Ala)4-X-(Ala)4 in solution, where X is the amino acid of interest. Further improvements include optimizing nonbonded interactions between charged residues to reproduce QM interaction energies of the charged-protein model compounds and experimental osmotic pressures. Validation of the optimized Drude protein FF includes MD simulations of a collection of peptides and proteins including β-sheet structures, as well as transmembrane ion channels. Results showed that the updated Drude-2019 protein FF yields smaller overall root-mean-square differences of proteins as compared to the additive CHARMM36m and Drude-2013 FFs as well as similar or improved agreement with experimental NMR properties, allowing for long time scale simulation studies of proteins and more complex biomolecular systems in conjunction with the remainder of the Drude polarizable FF.

Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids.

Computational simulations of ionic liquid solutions have become a useful tool to investigate various physical, chemical and catalytic properties of systems involving these solvents. Classical molecular dynamics and hybrid quantum mechanical/molecular mechanical (QM/MM) calculations of IL systems have provided significant insights at the atomic level. Here, we present a review of the development and application of the multipolar and polarizable force field AMOEBA for ionic liquid systems, termed AMOEBA–IL. The parametrization approach for AMOEBA–IL relies on the reproduction of total quantum mechanical (QM) intermolecular interaction energies and QM energy decomposition analysis. This approach has been used to develop parameters for imidazolium– and pyrrolidinium–based ILs coupled with various inorganic anions. AMOEBA–IL has been used to investigate and predict the properties of a variety of systems including neat ILs and IL mixtures, water exchange reactions on lanthanide ions in IL mixtures, IL–based liquid–liquid extraction, and effects of ILs on an aniline protection reaction.

The phenoxyl group-modulated interplay of cation-π and σ-type interactions in the alkali metal series.

The interaction potential energy surfaces (IPESs) of four alkaline metal cations (Na+, K+, Rb+ and Cs+) complexed with phenol and catechol were explored by accurate ab initio calculations to investigate the interplay of different noncovalent interactions and their behavior along the alkali metal series and upon -OH substitution. Selected one-dimensional interaction energy curves revealed two different minimum energy configurations for all phenol- and catechol-metal complexes, characterized either by cation-π or σ-type interactions. For each investigated complex several two-dimensional IPES maps were also computed, exploiting the computational advantages of the MP2mod approach. The size of the alkali cation was found to play a similar role in modulating both kinds of complexes, as the interaction strength always decreases along the metal series, from Na+ to Cs+. Conversely, the number of hydroxyl substituents markedly affected cation-π complexes vs. σ-type ones. As a most relevant finding, in catechol-metal complexes the strength of cation-π interactions is around half that of the σ-type ones. It is argued that the combined effect of cation dimensions and hydroxyl substitution in catechol-Na+ complexes makes σ-type configurations remarkably more stable and easily accessible than cation-π ones. Besides shedding new light on the origin of biological phenomena connected with underwater adhesion, the quantum mechanical interaction energy database provided herein may offer a useful reference for tuning accurate force fields, suitable for molecular dynamics simulations, where environmental effects might be also taken into account.

Improved Modeling of Cation-π and Anion-Ring Interactions Using the Drude Polarizable Empirical Force Field for Proteins.

Cation-π interactions are noncovalent interactions between a π-electron system and a positively charged ion that are regarded as a strong noncovalent interaction and are ubiquitous in biological systems. Similarly, though less studied, anion-ring interactions are present in proteins along with in-plane interactions of anions with aromatic rings. As these interactions are between a polarizing ion and a polarizable π system, the accuracy of the treatment of these interactions in molecular dynamics (MD) simulations using additive force fields (FFs) may be limited. In the present work, to allow for a better description of ion-π interactions in proteins in the Drude-2013 protein polarizable FF, we systematically optimized the parameters for these interactions targeting model compound quantum mechanical (QM) interaction energies with atom pair-specific Lennard-Jones parameters along with virtual particles as selected ring centroids introduced to target the QM interaction energies and geometries. Subsequently, MD simulations were performed on a series of protein structures where ion-π pairs occur to evaluate the optimized parameters in the context of the Drude-2013 FF. The resulting FF leads to a significant improvement in reproducing the ion-π pair distances observed in experimental protein structures, as well as a smaller root-mean-square differences and fluctuations of the overall protein structures from experimental structures. Accordingly, the optimized Drude-2013 protein polarizable FF is suggested for use in MD simulations of proteins where cation-π and anion-ring interactions are critical. © 2019 Wiley Periodicals, Inc.

Polarizable force field for RNA based on the classical drude oscillator.

RNA molecules are highly dynamic and capable of adopting a wide range of complex, folded structures. The factors driving the folding and dynamics of these structures are dependent on a balance of base pairing, hydration, base stacking, ion interactions, and the conformational sampling of the 2'-hydroxyl group in the ribose sugar. The representation of these features is a challenge for empirical force fields used in molecular dynamics simulations. Toward meeting this challenge, the inclusion of explicit electronic polarization is important in accurately modeling RNA structure. In this work, we present a polarizable force field for RNA based on the classical Drude oscillator model, which represents electronic degrees of freedom via negatively charged particles attached to their parent atoms by harmonic springs. Beginning with parametrization against quantum mechanical base stacking interaction energy and conformational energy data, we have extended the Drude-2017 nucleic acid force field to include RNA. The conformational sampling of a range of RNA sequences were used to validate the force field, including canonical A-form RNA duplexes, stem-loops, and complex tertiary folds that bind multiple Mg2+ ions. Overall, the Drude-2017 RNA force field reproduces important properties of these structures, including the conformational sampling of the 2'-hydroxyl and key interactions with Mg2+ ions. © 2018 Wiley Periodicals, Inc.

Improved Modeling of Halogenated Ligand-Protein Interactions Using the Drude Polarizable and CHARMM Additive Empirical Force Fields.

Halogenated ligands can participate in nonbonding interactions with proteins via halogen bond (XB) or halogen-hydrogen bond donor (X-HBD) interactions. In the context of molecular dynamics (MD) simulations, the accuracy of the simulations depends strongly on the force field (FF) used. To ensure good reproduction of XB and X-HBD interactions with proteins, we optimized the previously developed additive CHARMM36/CHARMM General force field (CGenFF) and Drude polarizable force field by including atom pair-specific Lennard-Jones parameters for aromatic halogen-protein interactions. The optimization targeted quantum mechanical interaction energy surfaces with the developed parameters then examined for their ability to reproduce experimental halogen-containing ligand-protein interactions in MD simulations. The calculated halogenated ligand interaction geometries were in good overall agreement with the experimental crystal data for both the polarizable and additive FFs, showing that these models can accurately treat both XB and X-HBD interactions. Analysis of the ligand-protein interactions shows significant contributions of polarizability to binding occurring in the Drude FF, with self-polarization energy making both favorable and unfavorable contributions to binding. Further analysis of the dipole moments from aqueous solution to protein indicates the polarizable FF accounts for subtle changes of the environment of the ligands that can impact binding. The present work demonstrates the utility of the updated additive CHARMM36/CGenFF and polarizable Drude FFs for the study of halogenated ligand-protein interactions in computer-aided drug design.

Application of a simple quantum chemical approach to ligand fragment scoring for Trypanosoma brucei pteridine reductase\xa01 inhibition

There is a need for improved and generally applicable scoring functions for fragment-based approaches to ligand design. Here, we evaluate the performance of a computationally efficient model for inhibitory activity estimation, which is composed only of multipole electrostatic energy and dispersion energy terms that approximate long-range ab initio quantum mechanical interaction energies. We find that computed energies correlate well with inhibitory activity for a compound series with varying substituents targeting two subpockets of the binding site of Trypanosoma brucei pteridine reductase 1. For one subpocket, we find that the model is more predictive for inhibitory activity than the ab initio interaction energy calculated at the MP2 level. Furthermore, the model is found to outperform a commonly used empirical scoring method. Finally, we show that the results for the two subpockets can be combined, which suggests that this simple nonempirical scoring function could be applied in fragment–based drug design.

Polarizable Force Field for DNA Based on the Classical Drude Oscillator: I. Refinement Using Quantum Mechanical Base Stacking and Conformational Energetics.

Empirical force fields seek to relate the configuration of a set of atoms to its energy, thus yielding the forces governing its dynamics, using classical physics rather than more expensive quantum mechanical calculations that are computationally intractable for large systems. Most force fields used to simulate biomolecular systems use fixed atomic partial charges, neglecting the influence of electronic polarization, instead making use of a mean-field approximation that may not be transferable across environments. Recent hardware and software developments make polarizable simulations feasible, and to this end, polarizable force fields represent the next generation of molecular dynamics simulation technology. In this work, we describe the refinement of a polarizable force field for DNA based on the classical Drude oscillator model by targeting quantum mechanical interaction energies and conformational energy profiles of model compounds necessary to build a complete DNA force field. The parametrization strategy employed in the present work seeks to correct weak base stacking in A- and B-DNA and the unwinding of Z-DNA observed in the previous version of the force field, called Drude-2013. Refinement of base nonbonded terms and reparametrization of dihedral terms in the glycosidic linkage, deoxyribofuranose rings, and important backbone torsions resulted in improved agreement with quantum mechanical potential energy surfaces. Notably, we expand on previous efforts by explicitly including Z-DNA conformational energetics in the refinement.

Quantum mechanical calculation of nanomaterial-ligand interaction energies by molecular fractionation with conjugated caps method

Molecular fractionation with conjugate caps (MFCC) method is introduced for the efficient estimation of quantum mechanical (QM) interaction energies between nanomaterial (carbon nanotube, fullerene, and graphene surface) and ligand (charged and neutral). In the calculations, nanomaterials are partitioned into small fragments and conjugated caps that are properly capped, and the interaction energies can be obtained through the summation of QM calculations of the fragments from which the contribution of the conjugated caps is removed. All the calculations were performed by density functional theory (DFT) and dispersion contributions for the attractive interactions were investigated by dispersion corrected DFT method. The predicted interaction energies by MFCC at each computational level are found to give excellent agreement with full system (FS) calculations with the mean energy deviation just a fractional kcal/mol. The accurate determination of nanomaterial-ligand interaction energies by MFCC suggests that it is an effective method for performing QM calculations on nanomaterial-ligand systems.

Full QM Calculation of RNA Energy Using Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps Method.

In this study, the electrostatically embedded generalized molecular fractionation with conjugate caps (concaps) method (EE-GMFCC) was employed for efficient linear-scaling quantum mechanical (QM) calculation of total energies of RNAs. In the EE-GMFCC approach, the total energy of RNA is calculated by taking a proper combination of the QM energy of each nucleotide-centric fragment with large caps or small caps (termed EE-GMFCC-LC and EE-GMFCC-SC, respectively) deducted by the energies of concaps. The two-body QM interaction energy between non-neighboring ribonucleotides which are spatially in close contact are also taken into account for the energy calculation. Numerical studies were carried out to calculate the total energies of a number of RNAs using the EE-GMFCC-LC and EE-GMFCC-SC methods at levels of the Hartree-Fock (HF) method, density functional theory (DFT), and second-order many-body perturbation theory (MP2), respectively. The results show that the efficiency of the EE-GMFCC-SC method is about 3 times faster than the EE-GMFCC-LC method with minimal accuracy sacrifice. The EE-GMFCC-SC method is also applied for relative energy calculations of 20 different conformers of two RNA systems using HF and DFT, respectively. Both single-point and relative energy calculations demonstrate that the EE-GMFCC method has deviations from the full system results of only a few kcal/mol.

Polarizable Multipole-Based Force Field for Aromatic Molecules and Nucleobases.

Aromatic molecules with π electrons are commonly involved in chemical and biological recognitions. For example, nucleobases play central roles in DNA/RNA structure and their interactions with proteins. The delocalization of the π electrons is responsible for the high polarizability of aromatic molecules. In this work, the AMOEBA force field has been developed and applied to 5 regular nucleobases and 12 aromatic molecules. The permanent electrostatic energy is expressed as atomic multipole interactions between atom pairs, and many-body polarization is accounted for by mutually induced atomic dipoles. We have systematically investigated aromatic ring stacking and aromatic-water interactions for nucleobases and aromatic molecules, as well as base–base hydrogen-bonding pair interactions, all at various distances and orientations. van der Waals parameters were determined by comparison to the quantum mechanical interaction energy of these dimers and fine-tuned using condensed phase simulation. By comparing to quantum mechanical calculations, we show that the resulting classical potential is able to accurately describe molecular polarizability, molecular vibrational frequency, and dimer interaction energy of these aromatic systems. Condensed phase properties, including hydration free energy, liquid density, and heat of vaporization, are also in good overall agreement with experimental values. The structures of benzene liquid phase and benzene-water solution were also investigated by simulation and compared with experimental and PDB structure derived statistical results.

Balancing the Interactions of Mg2+ in Aqueous Solution and with Nucleic Acid Moieties For a Polarizable Force Field Based on the Classical Drude Oscillator Model.

Mg2+ ions are important in biological systems, particularly in stabilizing compact RNA folds. Mg2+ is strongly polarizing, and representing its interactions in heterogeneous environments is a challenge for empirical force field development. To date, the most commonly used force fields in molecular dynamics simulations utilize a pairwise-additive approximation for electrostatic interactions, which cannot account for the significant polarization response in systems containing Mg2+. In the present work, we refine the interactions of Mg2+ with water, Cl- ions, and nucleic acid moieties using a polarizable force field based on the classical Drude oscillator model. By targeting gas-phase quantum mechanical interaction energies and geometries of hydrated complexes, as well as condensed-phase osmotic pressure calculations, we present a model for Mg2+ that yields quantitative agreement with experimental measurements of water dissociation free energy and osmotic pressure across a broad range of concentrations. Notable is the direct modeling of steric repulsion between the water Drude oscillators and Mg2+ to treat the Pauli exclusion effects associated with overlap of the electron clouds of water molecules in the first hydration shell around Mg2+. Combined with the refined interactions with nucleic acid moieties, the present model represents a significant advancement in simulating nucleic acid systems containing Mg2+.

Parametrization of halogen bonds in the CHARMM general force field: Improved treatment of ligand–protein interactions

A halogen bond is a highly directional, non-covalent interaction between a halogen atom and another electronegative atom. It arises due to the formation of a small region of positive electrostatic potential opposite the covalent bond to the halogen, called the ‘sigma hole.’ Empirical force fields in which the electrostatic interactions are represented by atom-centered point charges cannot capture this effect because halogen atoms usually carry a negative charge and therefore interact unfavorably with other electronegative atoms. A strategy to overcome this problem is to attach a positively charged virtual particle to the halogen. In this work, we extend the additive CHARMM General Force Field (CGenFF) to include such interactions in model systems of phenyl-X, with X being Cl, Br or I including di- and trihalogenated species. The charges, Lennard-Jones parameters, and halogen-virtual particle distances were optimized to reproduce the orientation dependence of quantum mechanical interaction energies with water, acetone, and N-methylacetamide as well as experimental pure liquid properties and relative hydration free energies with respect to benzene. The resulting parameters were validated in molecular dynamics simulations on small-molecule crystals and on solvated protein–ligand complexes containing halogenated compounds. The inclusion of positive virtual sites leads to better agreement across experimental observables, including preservation of ligand binding poses as a direct result of the improved representation of halogen bonding.

On the origin of donor O-H bond weakening in phenol-water complexes.

Matrix isolation infrared spectroscopy has been used to investigate intermolecular interactions in a series of binary O-H⋯O hydrogen bonded phenol-water complexes where water is the common acceptor. The interaction at the binding site has been tuned by incorporating multiple fluorine substitutions at different aromatic ring sites of the phenol moiety. The spectral effects for the aforesaid chemical changes are manifested in the infrared spectra of the complexes as systematic increase in spectral shift of the phenolic O-H stretching fundamental (ΔνO-H). While νO-H bands of the monomers of all the fluorophenols appear within a very narrow frequency range, the increase in ΔνO-H of the complexes from phenol to pentafluorophenol is very large, nearly 90%. The observed values of ΔνO-H do not show a linear correlation with the total binding energies (ΔEb) of the complexes, expected according to Badger-Bauer rule. However, in the same ΔνO-H vs ΔEb plot, nice linear correlations are revealed if the complexes of ortho-fluorophenols are treated separately from their meta/para-substituted analogues. The observations imply that in spite of having the same binding site (O-H⋯O) and the same chemical identities (phenolic), the complexes of ortho and non-ortho fluorophenols do not belong, from the viewpoint of detailed molecular interactions, to a homologous series. Linear correlations of ΔνO-H are, however, observed with respect to the electrostatic component of ΔEb as well as the quantum mechanical charge transfer interaction energy (ECT). From quantitative viewpoint, the latter correlation along with the associated electronic structure parameters appears more satisfactory. It has also been noted that the observed ΔνO-H values of the complexes display a linear relationship with the aqueous phase pKa values of the respective phenol derivatives.

Polarizable Multipole-Based Force Field for Dimethyl and Trimethyl Phosphate.

Phosphate groups are commonly observed in biomolecules such as nucleic acids and lipids. Due to their highly charged and polarizable nature, modeling these compounds with classical force fields is challenging. Using quantum mechanical studies and liquid-phase simulations, the AMOEBA force field for dimethyl phosphate (DMP) ion and trimethyl phosphate (TMP) has been developed. On the basis of ab initio calculations, it was found that ion binding and the solution environment significantly impact both the molecular geometry and the energy differences between conformations. Atomic multipole moments are derived from MP2/cc-pVQZ calculations of methyl phosphates at several conformations with their chemical environments taken into account. Many-body polarization is handled via a Thole-style induction model using distributed atomic polarizabilities. van der Waals parameters of phosphate and oxygen atoms are determined by fitting to the quantum mechanical interaction energy curves for water with DMP or TMP. Additional stretch-torsion and angle-torsion coupling terms were introduced in order to capture asymmetry in P–O bond lengths and angles due to the generalized anomeric effect. The resulting force field for DMP and TMP is able to accurately describe both the molecular structure and conformational energy surface, including bond and angle variations with conformation, as well as interaction of both species with water and metal ions. The force field was further validated for TMP in the condensed phase by computing hydration free energy, liquid density, and heat of vaporization. The polarization behavior between liquid TMP and TMP in water is drastically different.

Collecting high-order interactions in an effective pairwise intermolecular potential using the hydrated ion concept: the hydration of Cf³⁺.

This work proposes a new methodology to build interaction potentials between a highly charged metal cation and water molecules. These potentials, which can be used in classical computer simulations, have been fitted to reproduce quantum mechanical interaction energies (MP2 and BP86) for a wide range of [M(H2O)n](m+)(H2O)ℓ clusters (n going from 6 to 10 and ℓ from 0 to 18). A flexible and polarizable water shell model (Mobile Charge Density of Harmonic Oscillator) has been coupled to the cation-water potential. The simultaneous consideration of poly-hydrated clusters and the polarizability of the interacting particles allows the inclusion of the most important many-body effects in the new polarizable potential. Applications have been centered on the californium, Cf(III) the heaviest actinoid experimentally studied in solution. Two different strategies to select a set of about 2000 structures which are used for the potential building were checked. Monte Carlo simulations of Cf(III)+500 H2O for three of the intermolecular potentials predict an aquaion structure with coordination number close to 8 and average R(Cf-O) in the range 2.43-2.48 Å, whereas the fourth one is closer to 9 with R(Cf-O) = 2.54 Å. Simulated EXAFS spectra derived from the structural Monte Carlo distribution compares fairly well with the available experimental spectrum for the simulations bearing 8 water molecules. An angular distribution similar to that of a square antiprism is found for the octa-coordination.

Fragment quantum mechanical calculation of proteins and its applications.

Conspectus The desire to study molecular systems that are much larger than what the current state-of-the-art ab initio or density functional theory methods could handle has naturally led to the development of novel approximate methods, including semiempirical approaches, reduced-scaling methods, and fragmentation methods. The major computational limitation of ab initio methods is the scaling problem, because the cost of ab initio calculation scales nth power or worse with system size. In the past decade, the fragmentation approach based on chemical locality has opened a new door for developing linear-scaling quantum mechanical (QM) methods for large systems and for applications to large molecular systems such as biomolecules. The fragmentation approach is highly attractive from a computational standpoint. First, the ab initio calculation of individual fragments can be conducted almost independently, which makes it suitable for massively parallel computations. Second, the electron properties, such as density and energy, are typically combined in a linear fashion to reproduce those for the entire molecular system, which makes the overall computation scale linearly with the size of the system. In this Account, two fragmentation methods and their applications to macromolecules are described. They are the electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) method and the automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) approach. The EE-GMFCC method is developed from the MFCC approach, which was initially used to obtain accurate protein-ligand QM interaction energies. The main idea of the MFCC approach is that a pair of conjugate caps (concaps) is inserted at the location where the subsystem is divided by cutting the chemical bond. In addition, the pair of concaps is fused to form molecular species such that the overcounted effect from added concaps can be properly removed. By introducing the electrostatic embedding field in each fragment calculation and two-body interaction energy correction on top of the MFCC approach, the EE-GMFCC method is capable of accurately reproducing the QM molecular properties (such as the dipole moment, electron density, and electrostatic potential), the total energy, and the electrostatic solvation energy from full system calculations for proteins. On the other hand, the AF-QM/MM method was used for the efficient QM calculation of protein nuclear magnetic resonance (NMR) parameters, including the chemical shift, chemical shift anisotropy tensor, and spin-spin coupling constant. In the AF-QM/MM approach, each amino acid and all the residues in its vicinity are automatically assigned as the QM region through a distance cutoff for each residue-centric QM/MM calculation. Local chemical properties of the central residue can be obtained from individual QM/MM calculations. The AF-QM/MM approach precisely reproduces the NMR chemical shifts of proteins in the gas phase from full system QM calculations. Furthermore, via the incorporation of implicit and explicit solvent models, the protein NMR chemical shifts calculated by the AF-QM/MM method are in excellent agreement with experimental values. The applications of the AF-QM/MM method may also be extended to more general biological systems such as DNA/RNA and protein-ligand complexes.

Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps Method for Full Quantum Mechanical Calculation of Protein Energy

An electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) method is developed for efficient linear-scaling quantum mechanical (QM) calculation of protein energy. This approach is based on our previously proposed GMFCC/MM method (He; et al. J. Chem. Phys. 2006, 124, 184703), In this EE-GMFCC scheme, the total energy of protein is calculated by taking a linear combination of the QM energy of the neighboring residues and the two-body QM interaction energy between non-neighboring residues that are spatially in close contact. All the fragment calculations are embedded in a field of point charges representing the remaining protein environment, which is the major improvement over our previous GMFCC/MM approach. Numerical studies are carried out to calculate the total energies of 18 real three-dimensional proteins of up to 1142 atoms using the EE-GMFCC approach at the HF/6-31G* level. The overall mean unsigned error of EE-GMFCC for the 18 proteins is 2.39 kcal/mol with reference to the full system HF/6-31G* energies. The EE-GMFCC approach is also applied for proteins at the levels of the density functional theory (DFT) and second-order many-body perturbation theory (MP2), also showing only a few kcal/mol deviation from the corresponding full system result. The EE-GMFCC method is linear-scaling with a low prefactor, trivially parallel, and can be readily applied to routinely perform structural optimization of proteins and molecular dynamics simulation with high level ab initio electronic structure theories.