Enabling ab initio Hessian and frequency calculations of large molecules

Constructing Potential Energy Surface with Correlated Theory for Dipeptides Using Molecular Tailoring Approach

We demonstrate a cost-effective alternative employing the fragment-based molecular tailoring approach (MTA) for building the potential energy surface (PES) for two dipeptides viz. alanine-alanine and alanine-proline employing correlated theory, with augmented Dunning basis sets. About 1369 geometries are generated for each test dipeptide by systematically varying the dihedral angles Φ ${{\rm{\Phi }}}$ and Ψ ${{{\Psi }}}$ . These conformational geometries are partially optimized by relaxing all the other Z-matrix parameters, fixing the values of Φ ${{\rm{\Phi }}}$ and Ψ ${{{\Psi }}}$ . The MP2 level PES is constructed from the MTA-energies of chemically intact geometries using minimal hardware. The fidelity of MP2/aug-cc-pVDZ level PES is brought out by comparing it with its full calculation counterpart. Further, we bring out the power of the method by reporting the MTA-based CCSD/aug-cc-pVDZ level PES for these two dipeptides containing 498 and 562 basis functions respectively.

Molecular Tailoring Approach for the Estimation of Intramolecular Hydrogen Bond Energy

Hydrogen bonds (HBs) play a crucial role in many physicochemical and biological processes. Theoretical methods can reliably estimate the intermolecular HB energies. However, the methods for the quantification of intramolecular HB (IHB) energy available in the literature are mostly empirical or indirect and limited only to evaluating the energy of a single HB. During the past decade, the authors have developed a direct procedure for the IHB energy estimation based on the molecular tailoring approach (MTA), a fragmentation method. This MTA-based method can yield a reliable estimate of individual IHB energy in a system containing multiple H-bonds. After explaining and illustrating the methodology of MTA, we present its use for the IHB energy estimation in molecules and clusters. We also discuss the use of this method by other researchers as a standard, state-of-the-art method for estimating IHB energy as well as those of other noncovalent interactions.

Efficient Analytic Second Derivative of Electrostatic Embedding QM/MM Energy: Normal Mode Analysis of Plant Cryptochrome

Analytic second derivatives of electrostatic embedding (EE) quantum mechanics/molecular mechanics (QM/MM) energy are important for performing vibrational analysis and simulating vibrational spectra of quantum systems interacting with an environment represented as a classical electrostatic potential. The main bottleneck of EE-QM/MM second derivatives is the solution of coupled perturbed equations for each MM atom perturbation. Here, we exploit the Q-vector method [J. Chem. Phys., 2019, 151, 041102] to workaround this bottleneck. We derive the full analytic second derivative of the EE-QM/MM energy, which allows us to compute QM, MM, and QM-MM Hessian blocks in an efficient and easy to implement manner. To show the capabilities of our method, we compute the normal modes for the full Arabidopsis thaliana plant cryptochrome. We show that the flavin adenine dinucleotide vibrations (QM subsystem) strongly mix with protein modes. We compute approximate vibronic couplings for the lowest bright transition, from which we extract spectral densities and the homogeneous broadening of FAD absorption spectrum in protein using vibrationally resolved electronic spectrum simulations.

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.

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.

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.

Harmonizing accuracy and efficiency: A pragmatic approach to fragmentation of large molecules

Fragmentation methods offer an attractive alternative for ab initio treatment of large molecules and molecular clusters. However, balancing the accuracy and efficiency of these methods is a tight-rope-act. With this in view, we present an algorithm for automatic molecular fragmentation within Molecular Tailoring Approach (MTA) achieving this delicate balance. The automated code is tested out on a variety of molecules and clusters at the Hartree-Fock (HF)- and Møller-Plesset second order perturbation theory as well as density functional theory employing augmented Dunning basis sets. The results show remarkable accuracy and efficiency vis-à-vis the respective full calculations. Thus the present work forms an important step toward the development of an MTA-based black box code for implementation of HF as well as correlated quantum chemical calculations on large molecular systems.

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.

Analysis of Different Fragmentation Strategies on a Variety of Large Peptides: Implementation of a Low Level of Theory in Fragment-Based Methods Can Be a Crucial Factor

We have investigated the performance of two classes of fragmentation methods developed in our group (Molecules-in-Molecules (MIM) and Many-Overlapping-Body (MOB) expansion), to reproduce the unfragmented MP2 energies on a test set composed of 10 small to large biomolecules. They have also been assessed to recover the relative energies of different motifs of the acetyl(ala)18NH2 system. Performance of different bond-cutting environments and the use of Hartree-Fock and different density functionals (as a low level of theory) in conjunction with the fragmentation strategies have been analyzed. Our investigation shows that while a low level of theory (for recovering long-range interactions) may not be necessary for small peptides, it provides a very effective strategy to accurately reproduce the total and relative energies of larger peptides such as the different motifs of the acetyl(ala)18NH2 system. Employing M06-2X as the low level of theory, the calculated mean total energy deviation (maximum deviation) in the total MP2 energies for the 10 molecules in the test set at MIM(d=3.5Å), MIM(η=9), and MOB(d=5Å) are 1.16 (2.31), 0.72 (1.87), and 0.43 (2.02) kcal/mol, respectively. The excellent performance suggests that such fragment-based methods should be of general use for the computation of accurate energies of large biomolecular systems.

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.

Low energy isomers of (H2O)25 from a hierarchical method based on Monte Carlo temperature basin paving and molecular tailoring approaches benchmarked by MP2 calculations

We report new global minimum candidate structures for the (H2O)25 cluster that are lower in energy than the ones reported previously and correspond to hydrogen bonded networks with 42 hydrogen bonds and an interior, fully coordinated water molecule. These were obtained as a result of a hierarchical approach based on initial Monte Carlo Temperature Basin Paving sampling of the cluster's Potential Energy Surface with the Effective Fragment Potential, subsequent geometry optimization using the Molecular Tailoring Approach with the fragments treated at the second order Møller-Plesset (MP2) perturbation (MTA-MP2) and final refinement of the entire cluster at the MP2 level of theory. The MTA-MP2 optimized cluster geometries, constructed from the fragments, were found to be within <0.5 kcal/mol from the minimum geometries obtained from the MP2 optimization of the entire (H2O)25 cluster. In addition, the grafting of the MTA-MP2 energies yields electronic energies that are within <0.3 kcal/mol from the MP2 energies of the entire cluster while preserving their energy rank order. Finally, the MTA-MP2 approach was found to reproduce the MP2 harmonic vibrational frequencies, constructed from the fragments, quite accurately when compared to the MP2 ones of the entire cluster in both the HOH bending and the OH stretching regions of the spectra.

Analysis of hydrogen bond energies and hydrogen bonded networks in water clusters (H2O)20 and (H2O)25 using the charge-transfer and dispersion terms

The relationship of the charge-transfer and dispersion terms with the O–O length for every pair of hydrogen bonded water molecules in the isomers of (H₂O)₁₇–(H₂O)₂₁.

Appraisal of molecular tailoring approach for large clusters

High level ab initio investigations on molecular clusters are generally restricted to those of small size essentially due to the nonlinear scaling of corresponding computational cost. Molecular tailoring approach (MTA) is a fragmentation-based method, which offers an economical and efficient route for studying larger clusters. However, due to its approximate nature, the MTA-energies carry some errors vis-à-vis their full calculation counterparts. These errors in the MTA-energies are reduced by grafting the correction at a lower basis set (e.g., 6-31+G(d)) onto a higher basis set (e.g., aug-cc-pvdz or aug-cc-pvtz) calculation at MP2 level of theory. Further, better estimates of energies are obtained by making use of many-body interaction analysis. For this purpose, R-goodness (Rg) parameters for the three- and four-body interactions in a fragmentation scheme are proposed. The procedure employing grafting and many-body analysis has been tested out on molecular clusters of water, benzene, acetylene and carbon dioxide. It is found that for the fragmentation scheme having higher three- and four-body Rg-values, the errors in MTA-grafted energies are reduced typically to ~0.2 mH at MP2 level calculation. Coupled with the advantage in terms of computational resources and CPU time, the present method opens a possibility of accurate treatment of large molecular clusters.

Extrapolation to the Gold-Standard in Quantum Chemistry: Computationally Efficient and Accurate CCSD(T) Energies for Large Molecules Using an Automated Thermochemical Hierarchy

The CCSD(T) method is known as the gold-standard in quantum chemistry and has been the method of choice in quantum chemistry for over 20 years to obtain accurate bond energies and molecular properties. Its computational cost formally scales as the seventh power of the size of the system and can be prohibitive for large molecules. As part of our efforts to reduce the computational cost of the CCSD(T) method yet retain its accuracy, we present a simple, efficient, and user-friendly protocol to extrapolate to CCSD(T) energies in conjunction with MP2 energies. The method is based on the automated error-canceling thermochemical hierarchy previously developed by us called the Connectivity-Based Hierarchy (CBH). For a test set containing 30 diverse nonaromatic organic molecules and biomonomers, we obtain highly accurate extrapolated CCSD(T) energies (with a mean absolute error of only 0.2-0.3 kcal/mol with different basis-set). Additionally, the work also features the successful extrapolation to CCSD energies using a similar protocol.

Facilitating Minima Search for Large Water Clusters at the MP2 Level via Molecular Tailoring

Water clusters (H2O)20 and (H2O)25 are explored at the Møller-Plesset second-order perturbation (MP2) level of theory. Geometry optimization is carried out on favorable structures, initially generated by the temperature basin paving (TBP) method, utilizing the fragment-based molecular tailoring approach (MTA). MTA-based stabilization energies at the complete basis set limit are accurately estimated by grafting the energy correction using a smaller basis set. For prototypical cases, the minima are established via MTA-based vibrational frequency calculations at the MP2/aug-cc-pVDZ level. The potential of MTA in tackling large clusters is further demonstrated by performing geometry optimization at MP2/aug-cc-pVDZ starting with the global minimum of (H2O)30 reported by Monte Carlo (MC) and molecular dynamics (MD) investigations. The present study brings out the efficacy of MTA in performing computationally expensive ab initio calculations with minimal off-the-shelf hardware without significant loss of accuracy.