Many-body Expansion Research Articles

High-accuracy DMBE potential energy surface for CNO(A''4) and the rate coefficients for the C + NO reaction in the A'2, A''2, and A''4 states.

An accurate potential energy surface (PES) for the lowest lying A''4 state of the CNO system is presented based on explicitly correlated multi-reference configuration interaction calculations with quadruple zeta basis set (MRCI-F12/cc-pVQZ-F12). The abinitio energies are fitted using the double many-body expansion method, thus incorporating long-range energy terms that can accurately describe the electrostatic and dispersion interactions with physically motivated decaying functions. Together with the previously fitted lowest A'2 and A''2 states using the same theoretical framework, this constitutes a new set of PESs that are suitable to predict rate coefficients for all atom-diatom reactions of the CNO system. We use this set of PESs to calculate thermal rate coefficients for the C(P3) + NO(Π2) reaction and compare the temperature dependence and product branching ratios with experimental results. The comparison between theory and experiment is shown to be improved over previous theoretical studies. We highlight the importance of the long-range interactions for low-temperature rate coefficients.

Magnetic dipole operator from chiral effective field theory for many-body expansion methods

Magnetic dipole operator from chiral effective field theory for many-body expansion methods

Pitfalls in the n-mode representation of vibrational potentials.

Simulations of anharmonic vibrational motion rely on computationally expedient representations of the governing potential energy surface. The n-mode representation (n-MR)-effectively a many-body expansion in the space of molecular vibrations-is a general and efficient approach that is often used for this purpose in vibrational self-consistent field (VSCF) calculations and correlated analogues thereof. In the present analysis, a lack of convergence in many VSCF calculations is shown to originate from negative and unbound potentials at truncated orders of the n-MR expansion. For cases of strong anharmonic coupling between modes, the n-MR can both dip below the true global minimum of the potential surface and lead to effective single-mode potentials in VSCF that do not correspond to bound vibrational problems, even for bound total potentials. The present analysis serves mainly as a pathology report of this issue. Furthermore, this insight into the origin of VSCF non-convergence provides a simple, albeit ad hoc, route to correct the problem by "painting in" the full representation of groups of modes that exhibit these negative potentials at little additional computational cost. Somewhat surprisingly, this approach also reasonably approximates the results of the next-higher n-MR order and identifies groups of modes with particularly strong coupling. The method is shown to identify and correct problematic triples of modes-and restore SCF convergence-in two-mode representations of challenging test systems, including the water dimer and trimer, as well as protonated tropine.

Adaptive-Partitioning Multilayer Dynamics Simulations: 2. Implementations of the Permuted and Interpolated Adaptive-Partitioning Gradients.

Recently, an adaptive-partitioning multilayer Q1/Q2/MM method was proposed, where Q1 and Q2 denote, respectively, two distinct quantum-mechanical levels of theory and MM, the molecular-mechanical force fields. Such a multilayer model resembles the ONIOM (our own N-layered integrated molecular orbital and molecular mechanics) model by Morokuma and co-workers, but it is distinguished by on-the-fly reclassifying atoms to be Q1, Q2, or MM in dynamics simulations. To smoothly blend the levels of descriptions of the atoms, buffer zones are introduced between adjacent layers, and the energy is smoothly interpolated. In particular, the Q1/Q2 interaction energy was expressed in two different formalisms: permuted and interpolated adaptive-partitioning (PAP and IAP), respectively. While the PAP energy is based on a weighted many-body expansion, the IAP energy is derived via alchemical quantum calculations with interpolated Fock and overlap matrices. In this article, we examine in-depth the irregularities in the IAP energy near the boundary between the buffer and Q2 zones, which were found prominent in some calculations. These irregularities are due to basis-set linear dependencies, which can be effectively suppressed using a cutoff for the weighted atomic orbital coefficients. Furthermore, we derived and implemented the gradients for both PAP and IAP. Test calculations on a series of water cluster models show perfectly smooth gradients in PAP, while a minor discontinuity occurs in IAP gradients at the buffer/Q2 boundary. The energy and gradient discontinuities in IAP become smaller when moving the buffer/Q2 boundary further away from the Q1 center and when increasing the size of the basis sets used. Overall, those discontinuities are controllable, and possible ways to further diminish them are discussed.

Consistent density functional theory-based description of ion hydration through density-corrected many-body representations

Delocalization error constrains the accuracy of density functional theory in describing molecular interactions in ion–water systems. Using Na+ and Cl− in water as model systems, we calculate the effects of delocalization error in the SCAN functional for describing ion–water and water–water interactions in hydrated ions, and demonstrate that density-corrected SCAN (DC-SCAN) predicts n-body and interaction energies with an accuracy approaching coupled cluster theory. The performance of DC-SCAN is size-consistent, maintaining an accurate description of molecular interactions well beyond the first solvation shell. Molecular dynamics simulations at ambient conditions with many-body MB-SCAN(DC) potentials, derived from the many-body expansion, predict the solvation structure of Na+ and Cl− in quantitative agreement with reference data, while simultaneously reproducing the structure of liquid water. Beyond rationalizing the accuracy of density-corrected models of ion hydration, our findings suggest that our unified density-corrected MB formalism holds great promise for efficient DFT-based simulations of condensed-phase systems with chemical accuracy.

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 abinitio 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.

Coupled-Cluster Density-Based Many-Body Expansion.

While CCSD(T) is often considered the "gold standard" of computational chemistry, the scaling of its computational cost as N7 limits its applicability for large and complex molecular systems. In this work, we apply the density-based many-body expansion [ Int. J. Quantum Chem. 2020, 120, e26228] in combination with CCSD(T). The accuracy of this approach is assessed for neutral, protonated, and deprotonated water hexamers, as well as (H2O)16 and (H2O)17 clusters. For the neutral water clusters, we find that already with a density-based two-body expansion, we are able to approximate the supermolecular CCSD(T) energies within chemical accuracy (4 kJ/mol). This surpasses the accuracy that is achieved with a conventional, energy-based three-body expansion. We show that this accuracy can be maintained even when approximating the electron densities using Hartree-Fock instead of using coupled-cluster densities. The density-based many-body expansion thus offers a simple, resource-efficient, and highly parallelizable approach that makes CCSD(T)-quality calculations feasible where they would otherwise be prohibitively expensive.

Site-Specific Ionization Potentials and Electron Affinities in Large Molecular Systems at Coupled Cluster Level.

A many-body expansion of ionization potentials and electron affinities is developed based on a combination of the fragment molecular orbital method and equation-of-motion coupled-cluster (EOM-CC). In addition to site-specific values, obtained as one-body properties, pair and triple corrections are added to account for nonlocal EOM-CC contributions of the molecular environment of a chromophore. The developed method is applied to carboxylic acids, alkyl cations, a protein ubiquitin (Protein Data Bank ID 1UBQ), and a nano ribbon of white graphene elucidating the effect of environment on ionization potential and electron affinity.

The zero-point vibration contributions to energy and compressibility of light rare-gas solids

The zero-point vibration energy, Gruneisen coefficient, and zero-point vibration pressure of light rare-gas solids (LRGS)(He, Ne, Ar) are calculated within a certain range of atomic spacing using many-body expansion and ab initio calculation methods based on quantum and atomic cluster theories. The results are then compared with the previous theoretical calculations and experimental data. It is concluded that the smaller the molar volume of LRGS, the greater the crystal's zero-point vibration energy and vibration frequency. The zero-point vibration energy of LRGSis only 6% of the total many-body interaction energy but should not be neglected. The Grüneisen coefficient of LRGS increases as the molar volume increases. The zero-point vibration pressure is smaller, which is about 3% of the total pressure. The error caused by the pressure calculation is not significant, but when the molar volume is reduced to a fixed value, the zero-point vibration pressure must be considered and cannot be ignored. In addition, the calculation results of zero-point vibration pressure are compared to the previous theoretical calculations and experimental results.

Finite-temperature many-body perturbation theory for anharmonic vibrations: Recursions, algebraic reduction, second-quantized reduction, diagrammatic rules, linked-diagram theorem, finite-temperature self-consistent field, and general-order algorithm.

A unified theory is presented for finite-temperature many-body perturbation expansions of the anharmonic vibrational contributions to thermodynamic functions, i.e., the free energy, internal energy, and entropy. The theory is diagrammatically size-consistent at any order, as ensured by the linked-diagram theorem proved in this study, and, thus, applicable to molecular gases and solids on an equal footing. It is also a basis-set-free formalism, just like its underlying Bose-Einstein theory, capable of summing anharmonic effects over an infinite number of states analytically. It is formulated by the Rayleigh-Schrödinger-style recursions, generating sum-over-states formulas for the perturbation series, which unambiguously converges at the finite-temperature vibrational full-configuration-interaction limits. Two strategies are introduced to reduce these sum-over-states formulas into compact sum-over-modes analytical formulas. One is a purely algebraic method that factorizes each many-mode thermal average into a product of one-mode thermal averages, which are then evaluated by the thermal Born-Huang rules. Canonical forms of these rules are proposed, dramatically expediting the reduction process. The other is finite-temperature normal-ordered second quantization, which is fully developed in this study, including a proof of thermal Wick's theorem and the derivation of a normal-ordered vibrational Hamiltonian at finite temperature. The latter naturally defines a finite-temperature extension of size-extensive vibrational self-consistent field theory. These reduced formulas can be represented graphically as Feynman diagrams with resolvent lines, which include anomalous and renormalization diagrams. Two order-by-order and one general-order algorithms of computing these perturbation corrections are implemented and applied up to the eighth order. The results show no signs of Kohn-Luttinger-type nonconvergence.

Data efficiency and extrapolation trends in neural network interatomic potentials

Recently, key architectural advances have been proposed for neural network interatomic potentials (NNIPs), such as incorporating message-passing networks, equivariance, or many-body expansion terms. Although modern NNIP models exhibit small differences in test accuracy, this metric is still considered the main target when developing new NNIP architectures. In this work, we show how architectural and optimization choices influence the generalization of NNIPs, revealing trends in molecular dynamics (MD) stability, data efficiency, and loss landscapes. Using the 3BPA dataset, we uncover trends in NNIP errors and robustness to noise, showing these metrics are insufficient to predict MD stability in the high-accuracy regime. With a large-scale study on NequIP, MACE, and their optimizers, we show that our metric of loss entropy predicts out-of-distribution error and data efficiency despite being computed only on the training set. This work provides a deep learning justification for probing extrapolation and can inform the development of next-generation NNIPs.

SiS formation in the interstellar medium via SiH + S collisions

ABSTRACT One of the most important Si-bearing species in the intersellar medium is the SiS molecule. Thermal rate coefficients and other collisional properties are calculated for its formation via the title reaction using the quasi-classical trajectory method. An accurate representation of the HSiS potential energy surface is employed, which has been modelled from high-level ab initio calculations and a reliable description of long-range interactions as implied by the underlying double many-body expansion method. The calculated rate coefficients for the $\rm SiH + S \rightarrow SiS + H$ reaction can be modelled with k(T) = α(T/300)βe−γ/T where $\alpha =0.63\times 10^{-10}\, \rm cm^3\, s^{-1}$, β = -0.11, and $\gamma = 11.6\, \rm K$. This result is only slightly lower than that for SiS formation via Si + SH collisions. The contribution of each reaction mechanism and the rovibrational energy distributions of the nascent SiS molecule are also calculated. The title collision can also yield SH ($\rm SiH + S \rightarrow SH + Si$), but the corresponding rate coefficient is 20 to 27 times smaller than for SiS formation ($\alpha =0.025\times 10^{-10}\rm cm^3\, s^{-1}$, β =-0.13, and $\gamma = 9.38\, \rm K$). The role of intersections between excited electronic states is also discussed, based on novel calculations including eight electronic states.

An accurate and efficient fragmentation approach via the generalized many-body expansion for density matrices.

With relevant chemical space growing larger and larger by the day, the ability to extend computational tractability over that larger space is of paramount importance in virtually all fields of science. The solution we aim to provide here for this issue is in the form of the generalized many-body expansion for building density matrices (GMBE-DM) based on the set-theoretical derivation with overlapping fragments, through which the energy can be obtained by a single Fock build. In combination with the purification scheme and the truncation at the one-body level, the DM-based GMBE(1)-DM-P approach shows both highly accurate absolute and relative energies for medium-to-large size water clusters with about an order of magnitude better than the corresponding energy-based GMBE(1) scheme. Simultaneously, GMBE(1)-DM-P is about an order of magnitude faster than the previously proposed MBE-DM scheme [F. Ballesteros and K. U. Lao, J. Chem. Theory Comput. 18, 179 (2022)] and is even faster than a supersystem calculation without significant parallelization to rescue the fragmentation method. For even more challenging systems including ion-water and ion-pair clusters, GMBE(1)-DM-P also performs about 3 and 30 times better than the energy-based GMBE(1) approach, respectively. In addition, this work provides the first overlapping fragmentation algorithm with a robust and effective binning scheme implemented internally in a popular quantum chemistry software package. Thus, GMBE(1)-DM-P opens a new door to accurately and efficiently describe noncovalent clusters using quantum mechanics.

Orders-of-coupling representation achieved with a single neural network with optimal neuron activation functions and without nonlinear parameter optimization

Orders-of-coupling representations (representations of multivariate functions with low-dimensional functions that depend on subsets of original coordinates corresponding to different orders of coupling) are useful in many applications, for example, in computational chemistry and other applications, especially where integration is needed. Examples include N-mode approximations and many-body expansions. Such representations can be conveniently built with machine learning methods, and previously, methods building the lower-dimensional terms of such representations with neural networks [e.g. Comput. Phys. Commun. 180 (2009) 2002] and Gaussian process regressions [e.g. Mach. Learn. Sci. Technol. 3 (2022) 01LT02] were proposed. Here, we show that neural network models of orders-of-coupling representations can be easily built by using a recently proposed neural network with optimal neuron activation functions computed with a first-order additive Gaussian process regression [arXiv:2301.05567] and avoiding non-linear parameter optimization. Examples are given of representations of molecular potential energy surfaces.

MBX: A many-body energy and force calculator for data-driven many-body simulations.

Many-Body eXpansion (MBX) is a C++ library that implements many-body potential energy functions (PEFs) within the "many-body energy" (MB-nrg) formalism. MB-nrg PEFs integrate an underlying polarizable model with explicit machine-learned representations of many-body interactions to achieve chemical accuracy from the gas to the condensed phases. MBX can be employed either as a stand-alone package or as an energy/force engine that can be integrated with generic software for molecular dynamics and Monte Carlo simulations. MBX is parallelized internally using Open Multi-Processing and can utilize Message Passing Interface when available in interfaced molecular simulation software. MBX enables classical and quantum molecular simulations with MB-nrg PEFs, as well as hybrid simulations that combine conventional force fields and MB-nrg PEFs, for diverse systems ranging from small gas-phase clusters to aqueous solutions and molecular fluids to biomolecular systems and metal-organic frameworks.

Fixing the catastrophic breakdown of single reference coupled cluster theory for strongly correlated systems: Two paradigms toward the implicit inclusion of high-rank correlation with low-spin channels.

The dual exponential coupled cluster theory proposed by Tribedi et al.[J. Chem. Theory Comput. 16, 10, 6317-6328 (2020)] performs significantly better for a wide range of weakly correlated systems than the coupled cluster theory with singles and doubles excitations due to the implicit inclusion of high-rank excitations. The high-rank excitations are included through the action of a set of vacuum annihilating scattering operators that act non-trivially on certain correlated wavefunctions and are determined via a set of local denominators involving the energy difference between certain excited states. This often leads the theory to be prone to instabilities. In this paper, we show that restricting the correlated wavefunction, on which the scattering operators act, to be spanned by only the singlet-paired determinants can avoid catastrophic breakdown. For the first time, we present two nonequivalent approaches to arrive at the working equations, viz., the projective approach with sufficiency conditions and the amplitude form with many-body expansion. Although the effect of the triple excitation is quite small around molecular equilibrium geometry, this scheme leads to a better qualitative description of the energetics in the regions of strong correlation. With many pilot numerical applications, we have demonstrated the performance of the dual-exponential scheme with both the proposed solution strategies while restricting the excitation subspaces coupled to the corresponding lowest spin channels.

Benchmark coupled-cluster lattice energy of crystalline benzene and assessment of multi-level approximations in the many-body expansion.

The many-body expansion (MBE) is promising for the efficient, parallel computation of lattice energies in organic crystals. Very high accuracy should be achievable by employing coupled-cluster singles, doubles, and perturbative triples at the complete basis set limit [CCSD(T)/CBS] for the dimers, trimers, and potentially tetramers resulting from the MBE, but such a brute-force approach seems impractical for crystals of all but the smallest molecules. Here, we investigate hybrid or multi-level approaches that employ CCSD(T)/CBS only for the closest dimers and trimers and utilize much faster methods like Møller-Plesset perturbation theory (MP2) for more distant dimers and trimers. For trimers, MP2 is supplemented with the Axilrod-Teller-Muto (ATM) model of three-body dispersion. MP2(+ATM) is shown to be a very effective replacement for CCSD(T)/CBS for all but the closest dimers and trimers. A limited investigation of tetramers using CCSD(T)/CBS suggests that the four-body contribution is entirely negligible. The large set of CCSD(T)/CBS dimer and trimer data should be valuable in benchmarking approximate methods for molecular crystals and allows us to see that a literature estimate of the core-valence contribution of the closest dimers to the lattice energy using just MP2 was overbinding by 0.5 kJ mol-1, and an estimate of the three-body contribution from the closest trimers using the T0 approximation in local CCSD(T) was underbinding by 0.7 kJ mol-1. Our CCSD(T)/CBS best estimate of the 0K lattice energy is -54.01 kJ mol-1, compared to an estimated experimental value of -55.3 ± 2.2 kJ mol-1.

Classical Simulation of Short-Time Quantum Dynamics

Recent progress in the development of quantum technologies has enabled the direct investigation of the dynamics of increasingly complex quantum many-body systems. This motivates the study of the complexity of classical algorithms for this problem in order to benchmark quantum simulators and to delineate the regime of quantum advantage. Here, we present classical algorithms for approximating the dynamics of local observables and nonlocal quantities such as the Loschmidt echo, where the evolution is governed by a local Hamiltonian. For short times, their computational cost scales polynomially with the system size and the inverse of the approximation error. In the case of local observables, the proposed algorithm has a better dependence on the approximation error than algorithms based on the Lieb-Robinson bound. Our results use cluster-expansion techniques adapted to the dynamical setting, for which we give a novel proof of their convergence. This has important physical consequences besides our efficient algorithms. In particular, we establish a novel quantum speed limit, a bound on dynamical phase transitions, and a concentration bound for product states evolved for short times.Received 2 November 2022Accepted 25 April 2023DOI:https://doi.org/10.1103/PRXQuantum.4.020340Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.Published by the American Physical SocietyPhysics Subject Headings (PhySH)Research AreasQuantum information theoryQuantum statistical mechanicsPhysical SystemsLattice models in statistical physicsQuantum spin modelsTechniquesApproximation methods for many-body systemsComputational complexitySeries expansions & exact enumerationQuantum InformationStatistical Physics

Dynamics of C(3P) + OH(X 2Π) reaction on the new global HCO(X2A') potential energy surface.

A precise analytical potential energy surface (PES) of HCO(X2A') is fitted from a great quantity of abinitio energy points computed with the multi-reference configuration interaction method and aug-cc-pV(Q/5)Z basis sets. The whole energy points extrapolated to the complete basis set limit are fitted by the many-body expansion formula. The calculated topographic characteristics are analyzed and compared with the existing work to prove the precision of the present HCO(X2A') PES. By utilizing the time-dependent wave packet and quasi-classical trajectory methods, the reaction probabilities, integral cross sections, and rate constants are computed. The results are compared in detail with the former results carried out on the other PES. Moreover, the provided information on stereodynamics leads to an in-depth understanding of the role of collision energy in product distribution.

An accurate many-body expansion potential energy surface for HO2 (X 2A″) by extrapolation to the complete basis set limit and quantum dynamics of the related reaction O(3P) + OH(2Π)

A new global many-body expansion potential energy surface (PES) for the ground electronic state of HO2 is constructed. All ab initio energy points are employed using aug-cc-pV(Q+ d)Z and aug-cc-pV(5+ d)Z basis sets at the multi-reference conformational interaction level and then extrapolated to the complete basis set limit. The topographic characteristics of the global PES for HO2 (X 2A″) are discussed in detail and compared with experimental results, showing good agreement. We subsequently consider the O(3P) + OH(2Π) → O2(3 ) + H(2S) reaction using the time-dependent wave packet method and calculate dynamics behaviors, including the vibration–rotation distribution and reaction probabilities. The results indicate that the new PES is suitable for dynamics investigations of O(3P) + OH(2Π) → O2(3 ) + H(2S) and can be used as a component for PESs to construct larger oxygen/hydrogen containing systems.