Molecular forces, geometries, and frequencies by systematic molecular fragmentation including embedded charges

Simple Composite Approach to Efficiently Estimate Basis Set Limit CCSD(T) Harmonic Frequencies and Reaction Thermochemistry

We introduce a simple strategy that combines the G3(MP2) composite method and explicitly correlated coupled cluster CCSD(T)-F12 method to efficiently estimate complete basis set CCSD(T) molecular geometries and harmonic vibrational frequencies at the cost of a double-ζ basis set calculation. Based on a large test set of 61 neutral, ionic, and open-shell molecules, and additionally 31 molecules in the HFREQ2014 data set, we demonstrate that this composite strategy has an average accuracy of 2 cm-1 or better relative to complete basis set CCSD(T) values. Using this approach, we estimated 696 CCSD(T)/CBS reaction energies of small to medium-sized systems containing up to 6 heavy atoms and confirmed existing approximations that use small basis set density functional theory methods [e.g., M06-2X/6-31+G(d)] to calculate thermal contributions to reaction enthalpies and Gibbs free energies that are accurate to within 0.2 kcal mol-1 on average.

Seamless Multilayer─A Novel Total Energy Partition Scheme for Embedded and Hybrid Calculations

In this paper, we provide general formulation of a multilayer approach, covering both additive and subtractive quantum mechanics/molecular mechanics (QM/MM) as special cases. After that, we suggest a novel definition of QM/MM total energy based on the consideration of a system divided into three layers. In a simplified form, it is E = E Q M ( 1 + 2 ) - E Q M ( 2 ) + E M M ( 2 + 3 ) , where layers 1, 2, and 3 represent inner QM, outer QM, and classical MM regions, respectively. The novel formulation is also not limited by only QM/MM combination of methods─in fact, any computational methods can be combined in a hybrid calculation. In this paper, we call the new approach seamless multilayer. Test calculations performed for silica and boric oxide show that the new approach requires no QM/MM interface parameterization as well as no or very simple correction terms for boundary atoms. This can greatly facilitate QM/MM studies of covalent inorganic solids. However, test calculations of α-Al2O3 show that for ionic compounds, the new method requires some additional development.

Toward an extreme-scale electronic structure system

Electronic structure calculations have the potential to predict key matter transformations for applications of strategic technological importance, from drug discovery to material science and catalysis. However, a predictive physicochemical characterization of these processes often requires accurate quantum chemical modeling of complex molecular systems with hundreds to thousands of atoms. Due to the computationally demanding nature of electronic structure calculations and the complexity of modern high-performance computing hardware, quantum chemistry software has historically failed to operate at such large molecular scales with accuracy and speed that are useful in practice. In this paper, novel algorithms and software are presented that enable extreme-scale quantum chemistry capabilities with particular emphasis on exascale calculations. This includes the development and application of the multi-Graphics Processing Unit (GPU) library LibCChem 2.0 as part of the General Atomic and Molecular Electronic Structure System package and of the standalone Extreme-scale Electronic Structure System (EXESS), designed from the ground up for scaling on thousands of GPUs to perform high-performance accurate quantum chemistry calculations at unprecedented speed and molecular scales. Among various results, we report that the EXESS implementation enables Hartree-Fock/cc-pVDZ plus RI-MP2/cc-pVDZ/cc-pVDZ-RIFIT calculations on an ionic liquid system with 623 016 electrons and 146 592 atoms in less than 45 min using 27 600 GPUs on the Summit supercomputer with a 94.6% parallel efficiency.

Toward a general neural network force field for protein simulations: Refining the intramolecular interaction in protein

Molecular dynamics (MD) is an extremely powerful, highly effective, and widely used approach to understanding the nature of chemical processes in atomic details for proteins. The accuracy of results from MD simulations is highly dependent on force fields. Currently, molecular mechanical (MM) force fields are mainly utilized in MD simulations because of their low computational cost. Quantum mechanical (QM) calculation has high accuracy, but it is exceedingly time consuming for protein simulations. Machine learning (ML) provides the capability for generating accurate potential at the QM level without increasing much computational effort for specific systems that can be studied at the QM level. However, the construction of general machine learned force fields, needed for broad applications and large and complex systems, is still challenging. Here, general and transferable neural network (NN) force fields based on CHARMM force fields, named CHARMM-NN, are constructed for proteins by training NN models on 27 fragments partitioned from the residue-based systematic molecular fragmentation (rSMF) method. The NN for each fragment is based on atom types and uses new input features that are similar to MM inputs, including bonds, angles, dihedrals, and non-bonded terms, which enhance the compatibility of CHARMM-NN to MM MD and enable the implementation of CHARMM-NN force fields in different MD programs. While the main part of the energy of the protein is based on rSMF and NN, the nonbonded interactions between the fragments and with water are taken from the CHARMM force field through mechanical embedding. The validations of the method for dipeptides on geometric data, relative potential energies, and structural reorganization energies demonstrate that the CHARMM-NN local minima on the potential energy surface are very accurate approximations to QM, showing the success of CHARMM-NN for bonded interactions. However, the MD simulations on peptides and proteins indicate that more accurate methods to represent protein-water interactions in fragments and non-bonded interactions between fragments should be considered in the future improvement of CHARMM-NN, which can increase the accuracy of approximation beyond the current mechanical embedding QM/MM level.

Linear-Scaling Systematic Molecular Fragmentation Approach for Perturbation Theory and Coupled-Cluster Methods

The coupled-cluster (CC) singles and doubles with perturbative triples [CCSD(T)] method is frequently referred to as the “gold standard” of modern computational chemistry. However, the high computational cost of CCSD(T) [O(N⁷)], where N is the number of basis functions, limits its applications to small-sized chemical systems. To address this problem, efficient implementations of linear-scaling coupled-cluster methods, which employ the systematic molecular fragmentation (SMF) approach, are reported. In this study, we aim to do the following: (1) To achieve exact linear scaling and to obtain a pure ab initio approach, we revise the handling of nonbonded interactions in the SMF approach, denoted by LSSMF. (2) A new fragmentation algorithm, which yields smaller-sized fragments, that better fits high-level CC methods is introduced. (3) A modified nonbonded fragmentation scheme is proposed to enhance the existent algorithm. Performances of the LSSMF-CC approaches, such as LSSMF-CCSD(T), are compared with their canonical versions for a set of alkane molecules, C_nH_2n+2 (n = 6–10), which includes 142 molecules. Our results demonstrate that the LSSMF approach introduces negligible errors compared with the canonical methods; mean absolute errors (MAEs) are between 0.20 and 0.59 kcal mol^–1 for LSSMF(3,1)-CCSD(T). For a larger alkanes set (L12), C_nH_2n+2 (n = 50–70), the performance of LSSMF for the second-order perturbation theory (MP2) is investigated. For the L12 set, various bonded and nonbonded levels are considered. Our results demonstrate that the combination of bonded level 6 with nonbonded level 2, LSSMF(6,2), provides very accurate results for the MP2 method with a MAE value of 0.32 kcal mol^–1. The LSSMF(6,2) approach yields more than a 26-fold reduction in errors compared with LSSMF(3,1). Hence, we obtain substantial improvements over the original SMF approach. To illustrate the efficiency and applicability of the LSSMF-CCSD(T) approach, we consider an alkane molecule with 10,004 atoms. For this molecule, the LSSMF(3,1)-CCSD(T)/cc-pVTZ energy computation, on a Linux cluster with 100 nodes, 4 cores, and 5 GB of memory provided to each node, is performed just in ∼24 h. As a second test, we consider a biomolecular complex (PDB code: 1GLA), which includes 10,488 atoms, to assess the efficiency of the LSSMF approach. The LSSMF(3,1)-FNO–CCSD(T)/cc-pVTZ energy computation is completed in ∼7 days for the biomolecular complex. Hence, our results demonstrate that the LSSMF-CC approaches are very efficient. Overall, we conclude the following: (1) The LSSMF(m, n)-CCSD(T) methods can be reliably used for large-scale chemical systems, where the canonical methods are not computationally affordable. (2) The accuracy of bonded level 3 is not satisfactory for large chemical systems. (3) For high-accuracy studies, bonded level 5 (or higher) and nonbonded level 2 should be employed.

Effect of the Hydration Shell on the Carbonyl Vibration in the Ala-Leu-Ala-Leu Peptide

The vibrational spectrum of the Ala-Leu-Ala-Leu peptide in solution, computed from first-principles simulations, shows a prominent band in the amide I region that is assigned to stretching of carbonyl groups. Close inspection reveals combined but slightly different contributions by the three carbonyl groups of the peptide. The shift in their exact vibrational signature is in agreement with the different probabilities of these groups to form hydrogen bonds with the solvent. The central carbonyl group has a hydrogen bond probability intermediate to the other two groups due to interchanges between different hydrogen-bonded states. Analysis of the interaction energies of individual water molecules with that group shows that shifts in its frequency are directly related to the interactions with the water molecules in the first hydration shell. The interaction strength is well correlated with the hydrogen bond distance and hydrogen bond angle, though there is no perfect match, allowing geometrical criteria for hydrogen bonds to be used as long as the sampling is sufficient to consider averages. The hydrogen bond state of a carbonyl group can therefore serve as an indicator of the solvent's effect on the vibrational frequency.

Fantasy versus reality in fragment-based quantum chemistry

Since the introduction of the fragment molecular orbital method 20 years ago, fragment-based approaches have occupied a small but growing niche in quantum chemistry. These methods decompose a large molecular system into subsystems small enough to be amenable to electronic structure calculations, following which the subsystem information is reassembled in order to approximate an otherwise intractable supersystem calculation. Fragmentation sidesteps the steep rise (with respect to system size) in the cost of ab initio calculations, replacing it with a distributed cost across numerous computer processors. Such methods are attractive, in part, because they are easily parallelizable and therefore readily amenable to exascale computing. As such, there has been hope that distributed computing might offer the proverbial "free lunch" in quantum chemistry, with the entrée being high-level calculations on very large systems. While fragment-based quantum chemistry can count many success stories, there also exists a seedy underbelly of rarely acknowledged problems. As these methods begin to mature, it is time to have a serious conversation about what they can and cannot be expected to accomplish in the near future. Both successes and challenges are highlighted in this Perspective.

Accelerating the Calculation of Solute-Solvent Interaction Energies through Systematic Molecular Fragmentation

The method of systematic molecular fragmentation by annihilation (SMFA) is modified to apply to the interaction energy between a solute and solvent, where the solute is a pair of reacting molecules. For NH₃ + CH₃Cl as the solute, it is shown that SMFA can estimate (to chemical accuracy) the average binding energy of the solute in large water clusters containing up to 160 water molecules, at an appropriate level of electronic structure theory. The SMFA calculation can be carried out in a computation time that makes it feasible to estimate the solvation contribution to free energies of activation and reaction by ensemble averaging.

Variational Formulation of the Generalized Many-Body Expansion with Self-Consistent Charge Embedding: Simple and Correct Analytic Energy Gradient for Fragment-Based ab Initio Molecular Dynamics

The many-body expansion (MBE) and its extension to overlapping fragments, the generalized (G)MBE, constitute the theoretical basis for most fragment-based approaches for large-scale quantum chemistry. We reformulate the GMBE for use with embedding charges determined self-consistently from the fragment wave functions, in a manner that preserves the variational nature of the underlying self-consistent field method. As a result, the analytic gradient retains the simple "sum of fragment gradients" form that is often assumed in practice, sometimes incorrectly. This obviates (without approximation) the need to solve coupled-perturbed equations, and we demonstrate stable, fragment-based ab initio molecular dynamics simulations using this technique. Energy conservation fails when charge-response contributions to the Fock matrix are neglected, even while geometry optimizations and vibrational frequency calculations may yet be accurate. Stable simulations can be recovered by means of straightforward modifications introduced here, providing a general paradigm for fragment-based ab initio molecular dynamics.

Toward Building Protein Force Fields by Residue-Based Systematic Molecular Fragmentation and Neural Network

Accurate force fields are crucial for molecular dynamics investigation of complex biological systems. Building accurate protein force fields from quantum mechanical (QM) calculations is challenging due to the complexity of proteins and high computational costs of QM methods. In order to overcome these two difficulties, here we developed the residue-based systematic molecular fragmentation method to partition general proteins into only 20 types of amino acid dipeptides and one type of peptide bond at level 1. The total energy of proteins is the combination of the energies of these fragments. Each type of the fragments is then parametrized using neural network (NN) representation of the QM reference. Adopting NN representation can circumvent the limitation of the analytic form of classical molecular mechanics (MM) force fields. Using MM force fields as the baseline, our method adds NN representation of QM corrections at the length scale of amino acid dipeptides. We tested our force fields for both homogeneous and heterogeneous polypeptides. Energy and forces predicted by our force fields compare favorably with full QM calculations from tripeptides to decapeptides. Our development provides an efficient and accurate method of building protein force fields fully from ab initio QM calculations.

Simulations of infrared and Raman spectra in solution using the fragment molecular orbital method

Calculation of IR and Raman spectra in solution for large molecular systems made possible with analytic FMO/PCM Hessians.

Geometry optimizations with the incremental molecular fragmentation method

Nuclear energy gradients for the incremental molecular fragmentation (IMF) method presented in our previous work [Meitei OR, Heßelmann A, Molecular energies from an incremental fragmentation method, J Chem Phys 144(8):084109, 2016] have been derived. Using the second-order Møller–Plesset perturbation theory method to describe the bonded and nonbonded energy and gradient contributions and the uncorrelated Hartree–Fock method to describe the correction increment, it is shown that the IMF gradient can be easily computed by a sum of the underlying individual derivatives of the energy contributions. The performance of the method has been compared against the supermolecular method by optimizing the structures of a range of polyglycine molecules with up to 36 glycine residues in the chain. It is shown that with a sensible set of parameters used in the fragmentation the supermolecular structures can be fairly well reproduced. In a few cases the optimization with the IMF method leads to structures that differ from the supermolecular ones. It was found, however, that these are more stable geometries also on the supermolecular potential energy surface.

Force Field for Water Based on Neural Network

We developed a novel neural network-based force field for water based on training with high-level ab initio theory. The force field was built based on an electrostatically embedded many-body expansion method truncated at binary interactions. The many-body expansion method is a common strategy to partition the total Hamiltonian of large systems into a hierarchy of few-body terms. Neural networks were trained to represent electrostatically embedded one-body and two-body interactions, which require as input only one and two water molecule calculations at the level of ab initio electronic structure method CCSD/aug-cc-pVDZ embedded in the molecular mechanics water environment, making it efficient as a general force field construction approach. Structural and dynamic properties of liquid water calculated with our force field show good agreement with experimental results. We constructed two sets of neural network based force fields: nonpolarizable and polarizable force fields. Simulation results show that the nonpolarizable force field using fixed TIP3P charges has already behaved well, since polarization effects and many-body effects are implicitly included due to the electrostatic embedding scheme. Our results demonstrate that the electrostatically embedded many-body expansion combined with neural network provides a promising and systematic way to build next-generation force fields at high accuracy and low computational costs, especially for large systems.

Intramolecular interactions in sterically crowded hydrocarbon molecules

A molecular fragmentation method has been used to analyze the intramolecular interactions in the three molecules coupled diamantane, hexaphenylethane, and all-meta-tert-butyl substituted hexaphenylethane. The significance of these systems lies in the fact, that steric crowding effects enable a stabilization of the central carbon bond that possesses an extended length (1.6 to 1.7 Å) beyond conventional carbon-carbon bonds due to the steric repulsion of the attached hydrocarbon groups. The total stability of these molecules therefore depends on a delicate balance between attractive interaction forces on the one hand and on repulsive forces on the other hand. We have quantified the different interaction energy contributions using symmetry-adapted perturbation theory based on a density functional theory description of the monomers. It has been found that the attractive dispersion interactions increase more strongly with the level of crowding in the systems than the counteracting exchange interactions. This shows that steric crowding effects can have a significant impact on the structure and stability of large and branched molecules. © 2017 Wiley Periodicals, Inc.

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.

Vibrational infrared and Raman spectra of polypeptides: Fragments-in-fragments within molecular tailoring approach

The present work reports the calculation of vibrational infrared (IR) and Raman spectra of large molecular systems employing molecular tailoring approach (MTA). Further, it extends the grafting procedure for the accurate evaluation of IR and Raman spectra of large molecular systems, employing a new methodology termed as Fragments-in-Fragments (FIF), within MTA. Unlike the previous MTA-based studies, the accurate estimation of the requisite molecular properties is achieved without performing any full calculations (FC). The basic idea of the grafting procedure is implemented by invoking the nearly basis-set-independent nature of the MTA-based error vis-à-vis the respective FCs. FIF has been tested out for the estimation of the above molecular properties for three isomers, viz., β-strand, 310- and α-helix of acetyl(alanine)nNH2 (n = 10, 15) polypeptides, three conformers of doubly protonated gramicidin S decapeptide and trpzip2 protein (PDB id: 1LE1), respectively, employing BP86/TZVP, M06/6-311G**, and M05-2X/6-31G** levels of theory. For most of the cases, a maximum difference of 3 cm(-1) is achieved between the grafted-MTA frequencies and the corresponding FC values. Further, a comparison of the BP86/TZVP level IR and Raman spectra of α-helical (alanine)20 and its N-deuterated derivative shows an excellent agreement with the existing experimental spectra. In view of the requirement of only MTA-based calculations and the ability of FIF to work at any level of theory, the current methodology provides a cost-effective solution for obtaining accurate spectra of large molecular systems.

Microsolvation within the Systematic Molecular Fragmentation by Annihilation Approach

We have applied the systematic molecular fragmentation by annihilation (SMFA) fragmentation technique to glycine and DNA base pairs in water clusters, systems for which explicit solvation is believed to be important. The SMFA method was found to be capable of describing the structures, especially in handling the complexity of hydrogen bonding, with energies produced being comparable with those from full molecule results. Thus, the ability to break down large calculations into a manageable time without loss of accuracy shows promise for application to real biological systems for which these effects are relevant.

Molecular energies from an incremental fragmentation method

The systematic molecular fragmentation method by Collins and Deev [J. Chem. Phys. 125, 104104 (2006)] has been used to calculate total energies and relative conformational energies for a number of small and extended molecular systems. In contrast to the original approach by Collins, we have tested the accuracy of the fragmentation method by utilising an incremental scheme in which the energies at the lowest level of the fragmentation are calculated on an accurate quantum chemistry level while lower-cost methods are used to correct the low-level energies through a high-level fragmentation. In this work, the fragment energies at the lowest level of fragmentation were calculated using the random-phase approximation (RPA) and two recently developed extensions to the RPA while the incremental corrections at higher levels of the fragmentation were calculated using standard density functional theory (DFT) methods. The complete incremental fragmentation method has been shown to reproduce the supermolecule results with a very good accuracy, almost independent on the molecular type, size, or type of decomposition. The fragmentation method has also been used in conjunction with the DFT-SAPT (symmetry-adapted perturbation theory) method which enables a breakdown of the total nonbonding energy contributions into individual interaction energy terms. Finally, the potential problems of the method connected with the use of capping hydrogen atoms are analysed and two possible solutions are supplied.

Accurate vibrational spectra via molecular tailoring approach: a case study of water clusters at MP2 level

In spite of the recent advents in parallel algorithms and computer hardware, high-level calculation of vibrational spectra of large molecules is still an uphill task. To overcome this, significant effort has been devoted to the development of new algorithms based on fragmentation methods. The present work provides the details of an efficient and accurate procedure for computing the vibrational spectra of large clusters employing molecular tailoring approach (MTA). The errors in the Hessian matrix elements and dipole derivatives arising due to the approximation nature of MTA are reduced by grafting the corrections from a smaller basis set. The algorithm has been tested out for obtaining vibrational spectra of neutral and charged water clusters at Møller-Plesset second order level of theory, and benchmarking them against the respective full calculation (FC) and/or experimental results. For (H2O)16 clusters, the estimated vibrational frequencies are found to differ by a maximum of 2 cm(-1) with reference to the corresponding FC values. Unlike the FC, the MTA-based calculations including grafting procedure can be performed on a limited hardware, yet take a fraction of the FC time. The present methodology, thus, opens a possibility of the accurate estimation of the vibrational spectra of large molecular systems, which is otherwise impossible or formidable.

Pair-Pair Approximation to the Generalized Many-Body Expansion: An Alternative to the Four-Body Expansion for ab Initio Prediction of Protein Energetics via Molecular Fragmentation

We introduce a "pair-pair" approximation to the generalized many-body expansion (pp-GMBE) as an approximation to a traditional four-body expansion, the latter of which is accurate but quickly becomes numerically unstable and ultimately intractable as the number of "bodies" (fragments) increases. The pp-GMBE method achieves a good balance between accuracy and efficiency by defining significant fragment pairs and then fragment quartets. An efficient fragmentation scheme is introduced for proteins such that the largest subsystems contain about 60 atoms. Application of the pp-GMBE method to proteins with as many as 70 amino acids (1142 atoms) reveals that pp-GMBE energies are quite faithful to those obtained when the same level of density functional theory is applied to the entire macromolecule. When combined with embedding charges obtained from natural population analysis, the pp-GMBE approach affords absolute energies that differ by 1-3 kcal/mol from full supersystem results, but it yields conformational energy profiles that are practically indistinguishable from the supersystem calculation at the same level of theory.

Are fragment-based quantum chemistry methods applicable to medium-sized water clusters?

The GEBF method is demonstrated to be more accurate than the EE-MB method for medium-sized water clusters.

Fragment quantum chemical approach to geometry optimization and vibrational spectrum calculation of proteins

Geometry optimization and vibrational spectra (infrared and Raman spectra) calculations of proteins are carried out by a quantum chemical approach using the EE-GMFCC (electrostatically embedded generalized molecular fractionation with conjugate caps) method (J. Phys. Chem. A, 2013, 117, 7149).

H2 Adsorption in a Porous Crystal: Accurate First-Principles Quantum Simulation

A general method is presented for constructing, from ab initio quantum chemistry calculations, the potential energy surface (PES) for H2 absorbed in a porous crystalline material. The method is illustrated for the metal-organic framework material MOF-5. Rigid body quantum diffusion Monte Carlo simulations are used in the construction of the PES and to evaluate the quantum ground state of H2 in MOF-5, the zero-point energy, and the enthalpy of adsorption at 0 K.

Calculating nuclear magnetic resonance shieldings using systematic molecular fragmentation by annihilation

Systematic fragmentation accurately predicts theoretical chemical shieldings.