Convergence of the Many-Body Expansion for Energy and Forces for Classical Polarizable Models in the Condensed Phase

The Quality of Embedding Charges Is Critical for Convergence of Many-Body Expansions When BSSE Is Absent

There is conflicting evidence in the literature concerning the benefits of charge embedding on the convergence of many-body expansions (MBEs). Using a systematic series of water and ion-water clusters of varying size, this study indicates that the effects of charge embedding can be masked by basis set superposition error (BSSE). When BSSE is removed, this study demonstrates that charge embedding can significantly accelerate MBE convergence, where the electrostatically embedded two-body method, EE-MBE(2), can often yield accuracy close to the four-body method, MBE(4). Contrary to previous studies on smaller systems, this work shows that the performance of EE-MBE is highly sensitive to the charge model, with the best performance obtained when the natural population analysis (NPA) charge model is used and generated at the same level of theory used in the subsystem and supersystem calculations. It was demonstrated that the "3c" composite method, PBEh-3c, yields NPA atomic charges that are in excellent agreement with those obtained from supersystem density functional theory calculations. The linear-scaling X-Polarization method provides a more general approach to estimating these supersystem QM atomic charges, but its performance depends on how the fragments are defined.

Scalable generalized screening for high-order terms in the many-body expansion: Algorithm, open-source implementation, and demonstration

The many-body expansion lies at the heart of numerous fragment-based methods that are intended to sidestep the nonlinear scaling of ab initio quantum chemistry, making electronic structure calculations feasible in large systems. In principle, inclusion of higher-order n-body terms ought to improve the accuracy in a controllable way, but unfavorable combinatorics often defeats this in practice and applications with n ≥ 4 are rare. Here, we outline an algorithm to overcome this combinatorial bottleneck, based on a bottom-up approach to energy-based screening. This is implemented within a new open-source software application ("Fragme∩t"), which is integrated with a lightweight semi-empirical method that is used to cull subsystems, attenuating the combinatorial growth of higher-order terms in the graph that is used to manage the calculations. This facilitates applications of unprecedented size, and we report four-body calculations in (H2O)64 clusters that afford relative energies within 0.1 kcal/mol/monomer of the supersystem result using less than 10% of the unique subsystems. We also report n-body calculations in (H2O)20 clusters up to n = 8, at which point the expansion terminates naturally due to screening. These are the largest n-body calculations reported to date using ab initio electronic structure theory, and they confirm that high-order n-body terms are mostly artifacts of basis-set superposition error.

Many-Body Effects in Aqueous Systems: Synergies Between Interaction Analysis Techniques and Force Field Development

Interaction analysis techniques, including the many-body expansion (MBE), symmetry-adapted perturbation theory, and energy decomposition analysis, allow for an intuitive understanding of complex molecular interactions. We review these methods by first providing a historical context for the study of many-body interactions and discussing how nonadditivities emerge from Hamiltonians containing strictly pairwise-additive interactions. We then elaborate on the synergy between these interaction analysis techniques and the development of advanced force fields aimed at accurately reproducing the Born-Oppenheimer potential energy surface. In particular, we focus on ab initio-based force fields that aim to explicitly reproduce many-body terms and are fitted to high-level electronic structure results. These force fields generally incorporate many-body effects through (a) parameterization of distributed multipoles, (b) explicit fitting of the MBE, (c) inclusion of many-atom features in a neural network, and (d) coarse-graining of many-body terms into an effective two-body term. We also discuss the emerging use of the MBE to improve the accuracy and speed of ab initio molecular dynamics.

An extensive assessment of the performance of pairwise and many-body interaction potentials in reproducing ab initio benchmark binding energies for water clusters n = 2-25

We assess the performance of 7 pairwise additive (TIP3P, TIP4P, TIP4P-ice, TIP5P, OPC, SPC, SPC/E) and 8 families of many-body potentials (q-AQUA, HIPPO, AMOEBA, EFP, TTM, WHBB, MB-pol, MB-UCB) in reproducing high-level ab initio benchmark values, CCSD(T) or MP2 at the complete basis set (CBS) limit for the binding energy and the many-body expansion (MBE) of water clusters n = 2-11, 16-17, 20, 25. By including a large range of cluster sizes having dissimilar hydrogen bonding networks, we obtain an understanding of how these potentials perform for different hydrogen bonding arrangements that are mostly outside of their parameterization range. While it is appropriate to compare the results of ab initio based many-body potentials directly to the electronic binding energies (D_e's), the pairwise additive ones are compared to the enthalpies at T = 298 K, ΔH(298 K), as the latter class of force fields are parametrized to reproduce enthalpies (implicitly accounting for zero-point energy corrections) rather than binding energies. We find that all pairwise additive potentials considered overestimate the reference ΔH values for the n = 2-25 clusters by >13%. For the water dimer (n = 2) in particular, the errors are in the range 83-119% for the pairwise additive potentials studied since these are based on an effective rather than the true 2-body interaction specifically designed as a means of partially accounting for the missing many-body terms. This stronger 2-body interaction is achieved by an enhanced monomer dipole moment that mimics its increase from the gas phase monomer to the condensed phase value. Indeed, for cluster sizes n ≥ 4 the percent deviations become slightly smaller (albeit all exceeding 13%). In contrast, we find that the many-body potentials perform more accurately in reproducing the electronic binding energies (D_e's) throughout the entire cluster range (n = 2-25), all reproducing the ab initio benchmark binding energies within ±7% of the respective CBS values. We further assess the ability of a subset of the many-body potentials (MB-UCB, q-AQUA, MB-pol, and TTM2.1-F) to also reproduce the magnitude of the ab initio many-body energy terms for water cluster sizes n = 7, 10, 16 and 17. The potentials show an overall good agreement with the available benchmark values. However, we identify characteristic differences upon comparing the many-body terms at both the ab initio-optimized geometries and the respective potential-optimized geometries to the reference ab initio values. Additionally, by applying this analysis to a wide range of cluster sizes, trends in the MBE of the potentials with increasing cluster size can be identified. Finally, in an attempt to draw a parallel between the pairwise additive and many-body potentials, we report the analysis of the individual molecular dipole moments for water clusters with 1 to ∼4 solvation shells with the TTM2.1-F potential. We find that the internally solvated water molecules have in general a larger molecular dipole moment ranging from 2.6-3.0 D. This justifies the use of an enhanced, with respect to the gas-phase value, molecular dipole moment for the pairwise additive potentials, which is intended to fold in the many body terms into an effective (enhanced) pairwise interaction through the choice of the charges. These results have important implications for the development of future generations of efficient, transferable, and highly accurate classical interaction potentials in both the pairwise additive and many-body categories.

A Formulation of the Many-Body Expansion (MBE) for Periodic Systems: Application to Several Ice Phases

We introduce a new formulation of the many-body expansion (MBE) for periodic systems and apply it to 7 ice polymorphs (Ih, II, VIII, IX, XIII, XIV, and XV). This new formulation is built via a hierarchical procedure that connects gas-phase clusters that mimic unit cells over finite supercells to infinite solids. For periodic systems, the method is validated by showing that the lattice energies computed up to the 4-body in the MBE reproduce the lattice energies obtained using periodic boundary conditions with an Ewald summation for the 7 ice polymorphs. This development makes it possible to quantify, for the first time, the many-body contributions to the lattice energy of various ice polymorphs, which vary significantly among the 7 ice phases, amounting to between 7 and 24% of the total lattice energies. This development opens the door for obtaining insights into solid-state properties, while leveraging the computational benefits of the MBE.

Evaluation of molecular photophysical and photochemical properties using linear response time-dependent density functional theory with classical embedding: Successes and challenges

Time-dependent density functional theory (TDDFT) based approaches have been developed in recent years to model the excited-state properties and transition processes of the molecules in the gas-phase and in a condensed medium, such as in a solution and protein microenvironment or near semiconductor and metal surfaces. In the latter case, usually, classical embedding models have been adopted to account for the molecular environmental effects, leading to the multi-scale approaches of TDDFT/polarizable continuum model (PCM) and TDDFT/molecular mechanics (MM), where a molecular system of interest is designated as the quantum mechanical region and treated with TDDFT, while the environment is usually described using either a PCM or (non-polarizable or polarizable) MM force fields. In this Perspective, we briefly review these TDDFT-related multi-scale models with a specific emphasis on the implementation of analytical energy derivatives, such as the energy gradient and Hessian, the nonadiabatic coupling, the spin-orbit coupling, and the transition dipole moment as well as their nuclear derivatives for various radiative and radiativeless transition processes among electronic states. Three variations of the TDDFT method, the Tamm-Dancoff approximation to TDDFT, spin-flip DFT, and spin-adiabatic TDDFT, are discussed. Moreover, using a model system (pyridine-Ag20 complex), we emphasize that caution is needed to properly account for system-environment interactions within the TDDFT/MM models. Specifically, one should appropriately damp the electrostatic embedding potential from MM atoms and carefully tune the van der Waals interaction potential between the system and the environment. We also highlight the lack of proper treatment of charge transfer between the quantum mechanics and MM regions as well as the need for accelerated TDDFT modelings and interpretability, which calls for new method developments.

Molecular Dynamics Driven by the Many-Body Expansion (MBE-MD)

We present a protocol for classical and nuclear quantum dynamics, in which the energies and forces are generated by the many-body expansion (MBE), and apply it to water clusters using the TTM2.1-F and MB-Pol interaction potentials at various temperatures. We carry out MBE-molecular dynamics (MD) classical and nuclear quantum dynamical simulations, in which the energies and forces of the full system are approximated by the two-, three-, and four-body terms of the MBE, and compare the average potential and the vibrational density of states with the full simulation, i.e., the one for which no MBE is used. Our results indicate that the thermally averaged potential energy from the MBE up to the four-body term converges with near-identical behavior to the one from the full simulation. The three-body makes a substantial contribution (∼20%) to the energy, whereas the four-body is necessary for obtaining quantitatively accurate energetics and forces, albeit making a small contribution to each (∼2%). We further show that the harmonic frequencies are reproduced to within a few wavenumbers (cm^-1) at the four-body level and that the slowest modes to converge with the MBE rank are those involving the strongest hydrogen bonds. Anharmonicity exacerbates this effect, so that a four-body description of the energies and forces is needed to achieve accurate anharmonic vibrational frequencies in the hydrogen-bonded OH-stretching region. We also discuss the asymptotic scaling of the MBE-MD protocol with respect to the cost of the underlying potential energy evaluation, suggesting that electronic structure methods that scale at least as N⁴, N being the size of the system, are needed to result in savings over the traditional full MD simulation. We anticipate that the MBE-MD protocol can evolve into a powerful and practical method, which will allow for highly accurate ab initio MD simulations on a much broader range of molecular systems than can be currently handled.

Cluster many-body expansion: A many-body expansion of the electron correlation energy about a cluster mean field reference

The many-body expansion (MBE) is an efficient tool that has a long history of use for calculating interaction energies, binding energies, lattice energies, and so on. In the past, applications of MBE to correlation energy have been unfeasible for large systems, but recent improvements to computing resources have sparked renewed interest in capturing the correlation energy using the generalized nth order Bethe-Goldstone equation. In this work, we extend this approach, originally proposed for a Slater determinant, to a tensor product state (TPS) based wavefunction. By partitioning the active space into smaller orbital clusters, our approach starts from a cluster mean field reference TPS configuration and includes the correlation contribution of the excited TPSs using the MBE. This method, named cluster MBE (cMBE), improves the convergence of MBE at lower orders compared to directly doing a block-based MBE from a RHF reference. We present numerical results for strongly correlated systems, such as the one- and two-dimensional Hubbard models and the chromium dimer. The performance of the cMBE method is also tested by partitioning the extended π space of several large π-conjugated systems, including a graphene nano-sheet with a very large active space of 114 electrons in 114 orbitals, which would require 10⁶⁶ determinants for the exact FCI solution.

From Intermolecular Interaction Energies and Observable Shifts to Component Contributions and Back Again: A Tale of Variational Energy Decomposition Analysis

Quantum chemistry in the form of density functional theory (DFT) calculations is a powerful numerical experiment for predicting intermolecular interaction energies. However, no chemical insight is gained in this way beyond predictions of observables. Energy decomposition analysis (EDA) can quantitatively bridge this gap by providing values for the chemical drivers of the interactions, such as permanent electrostatics, Pauli repulsion, dispersion, and charge transfer. These energetic contributions are identified by performing DFT calculations with constraints that disable components of the interaction. This review describes the second-generation version of the absolutely localized molecular orbital EDA (ALMO-EDA-II). The effects of different physical contributions on changes in observables such as structure and vibrational frequencies upon complex formation are characterized via the adiabatic EDA. Example applications include red- versus blue-shifting hydrogen bonds; the bonding and frequency shifts of CO, N₂, and BF bound to a [Ru(II)(NH₃)₅]^{2 +} moiety; and the nature of the strongly bound complexes between pyridine and the benzene and naphthalene radical cations. Additionally, the use of ALMO-EDA-II to benchmark and guide the development of advanced force fields for molecular simulation is illustrated with the recent, very promising, MB-UCB potential.

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.

Massively Parallel Implementation of Divide-and-Conquer Jacobi Iterations Using Particle-Mesh Ewald for Force Field Polarization

To accelerate the evaluation of the self-consistent polarization in large condensed-phase systems with polarizable force fields, the new divide-and-conquer Jacobi iterations (DC-JI) solver is adapted for periodic boundary conditions with particle-mesh Ewald and implemented in a massively parallel fashion within the Tinker-HP software package. DC-JI captures the mutual polarization of close-range interactions within subclusters of atoms using Cholesky decomposition and couples in the polarization effects between these clusters iteratively. Iterative convergence is accelerated with direct inversion of the iterative subspace (DIIS) extrapolation. Compared to widely used preconditioned conjugate gradient (PCG) or conventional Jacobi iterations (JI/DIIS) algorithms, DC-JI/DIIS solves the polarization equations ∼20-30% faster in protein systems ranging from ∼10,000-175,000 atoms run on hundreds of processor cores. This translates to ∼10-15% speed-ups in the number of nanoseconds of simulation time that can be achieved per day. Not only is DC-JI/DIIS faster than PCG, but it also gives more energetically robust solutions for a given convergence threshold. These improvements make numerically robust polarizable force field simulations more computationally tractable for chemical systems of interest.

Development of a Flexible-Monomer Two-Body Carbon Dioxide Potential and Its Application to Clusters up to (CO2 )13

A flexible-monomer two-body potential energy function was developed that approaches the high level CCSD(T)/CBS potential energy surface (PES) of carbon dioxide (CO₂ ) systems. This function was generated by fitting the electronic energies of unique CO₂ monomers and dimers to permutationally invariant polynomials. More than 200,000 CO₂ configurations were used to train the potential function. Comparisons of the PESs of six orientations of flexible CO₂ dimers were evaluated to demonstrate the accuracy of the potential. Furthermore, the potential function was used to determine the minimum energy structures of CO₂ clusters containing as many as 13 molecules. For isomers of (CO₂ )₃ , the potential demonstrated energetic agreement with the M06-2X functional and structural agreement of the B2PLYP-D functional at substantially reduced computational costs. A separate function, fit to MP2/aug-cc-pVDZ reference energies, was developed to directly compare the two-body potential to the ab initio MP2 level of theory. © 2017 Wiley Periodicals, Inc.

Accurate Classical Polarization Solution with No Self-Consistent Field Iterations

We present a new solution for classical polarization that does not require any self-consistent field iterations, the aspect of classical polarization that makes it computationally expensive. The new method builds upon our iEL/SCF Lagrangian scheme that defines a set of auxiliary induced dipoles whose original purpose was to serve as a time-reversible initial guess to the SCF solution of the set of real induced dipoles. In the new iEL/0-SCF approach the auxiliary dipoles now drive the time evolution of the real induced dipoles such that they stay close to the Born-Oppenheimer surface in order to achieve a truly SCF-less method. We show that the iEL/0-SCF exhibits no loss of simulation accuracy when analyzed across bulk water, low to high concentration salt solutions, and small solutes to large proteins in water. In addition, iEL/0-SCF offers significant computational savings over more expensive SCF calculations based on traditional 1 fs time step integration using symplectic integrators and is as fast as reversible reference system propagator algorithms with an outer 2 fs time step.

Advanced Potential Energy Surfaces for Molecular Simulation

Advanced potential energy surfaces are defined as theoretical models that explicitly include many-body effects that transcend the standard fixed-charge, pairwise-additive paradigm typically used in molecular simulation. However, several factors relating to their software implementation have precluded their widespread use in condensed-phase simulations: the computational cost of the theoretical models, a paucity of approximate models and algorithmic improvements that can ameliorate their cost, underdeveloped interfaces and limited dissemination in computational code bases that are widely used in the computational chemistry community, and software implementations that have not kept pace with modern high-performance computing (HPC) architectures, such as multicore CPUs and modern graphics processing units (GPUs). In this Feature Article we review recent progress made in these areas, including well-defined polarization approximations and new multipole electrostatic formulations, novel methods for solving the mutual polarization equations and increasing the MD time step, combining linear-scaling electronic structure methods with new QM/MM methods that account for mutual polarization between the two regions, and the greatly improved software deployment of these models and methods onto GPU and CPU hardware platforms. We have now approached an era where multipole-based polarizable force fields can be routinely used to obtain computational results comparable to state-of-the-art density functional theory while reaching sampling statistics that are acceptable when compared to that obtained from simpler fixed partial charge force fields.