Tensor Hypercontraction Second-Order Møller–Plesset Perturbation Theory: Grid Optimization and Reaction Energies

Tensor hypercontraction for fully self-consistent imaginary-time GF2 and GWSOX methods: Theory, implementation, and role of the Green's function second-order exchange for intermolecular interactions

We present an efficient MPI-parallel algorithm and its implementation for evaluating the self-consistent correlated second-order exchange term (SOX), which is employed as a correction to the fully self-consistent GW scheme called scGWSOX (GW plus the SOX term iterated to achieve full Green's function self-consistency). Due to the application of the tensor hypercontraction (THC) in our computational procedure, the scaling of the evaluation of scGWSOX is reduced from O(nτnAO5) to O(nτN2nAO2). This fully MPI-parallel and THC-adapted approach enabled us to conduct the largest fully self-consistent scGWSOX calculations with over 1100 atomic orbitals with only negligible errors attributed to THC fitting. Utilizing our THC implementation for scGW, scGF2, and scGWSOX, we evaluated energies of intermolecular interactions. This approach allowed us to circumvent issues related to reference dependence and ambiguity in energy evaluation, which are common challenges in non-self-consistent calculations. We demonstrate that scGW exhibits a slight overbinding tendency for large systems, contrary to the underbinding observed with non-self-consistent RPA. Conversely, scGWSOX exhibits a slight underbinding tendency for such systems. This behavior is both physical and systematic and is caused by exclusion-principle violating diagrams or corresponding corrections. Our analysis elucidates the role played by these different diagrams, which is crucial for the construction of rigorous, accurate, and systematic methods. Finally, we explicitly show that all perturbative fully self-consistent Green's function methods are size-extensive and size-consistent.

Tensor Hypercontraction of Cluster Perturbation Theory: Quartic Scaling Perturbation Series for the Coupled Cluster Singles and Doubles Ground-State Energies

Even though cluster perturbation theory has been shown to be a robust noniterative alternative to coupled cluster theory, it is still plagued by high order polynomial computational scaling and the storage of higher order tensors. We present a proof-of-concept strategy for implementing a cluster perturbation theory ground-state energy series for the coupled cluster singles and doubles energy with N4 computational scaling using tensor hypercontraction (THC). The reduction in computational scaling by two orders is achieved by decomposing two electron repulsion integrals, doubles amplitudes and multipliers, as well as selected double intermediates to the THC format. Using the outlined strategy, we showcase that the THC pilot implementations retain numerical accuracy to within 1 kcal/mol relative to corresponding conventional and density fitting implementations, and we empirically verify the N4 scaling.

Robust Tensor Hypercontraction of the Particle-Particle Ladder Term in Equation-of-Motion Coupled Cluster Theory

One method of representing a high-rank tensor as a (hyper-)product of lower-rank tensors is the tensor hypercontraction (THC) method of Hohenstein et al. This strategy has been found to be useful for reducing the polynomial scaling of coupled-cluster methods by representation of a four-dimensional tensor of electron-repulsion integrals in terms of five two-dimensional matrices. Pierce et al. have already shown that the application of a robust form of THC to the particle-particle ladder (PPL) term reduces the cost of this term in couple-cluster singles and doubles (CCSD) from O ( N 6 ) to O ( N 5 ) with negligible errors in energy with respect to the density-fitted variant. In this work, we have implemented the least-squares variant of THC (LS-THC) which does not require a nonlinear tensor factorization, including the robust form (R-LS-THC), for the calculation of the excitation and electron attachment energies using equation-of-motion coupled cluster methods EOMEE-CCSD and EOMEA-CCSD, respectively. We have benchmarked the effect of the R-LS-THC-PPL approximation on excitation energies using the comprehensive QUEST database and the accuracy of electron attachment energies using the NAB22 database. We find that errors on the order of 1 meV are achievable with a reduction in total calculation time of approximately 5 ×.

Local Second-Order Møller-Plesset Theory with a Single Threshold Using Orthogonal Virtual Orbitals: Theory, Implementation, and Assessment

It has long been clear that electron correlation methods exhibit unphysical compute scalings with molecular size, which has motivated the development of local correlation methods to discard effectively zero contributions in a controlled way to yield an approximate correlation energy. The ideal local correlation method should have a single numerical threshold that controls the dropping of terms with the ability to have that threshold set small enough so that the correlation energy is reproduced to enough significant figures such that the result is chemically identical. This work reports such a method for the second-order Møller-Plesset (MP2) theory. The theory, implementation, and testing of this local MP2 theory are reported. Thresholds ranging from 10-5 to 10-8 and basis sets ranging from split valence plus polarization through to quadruple-ζ are assessed for local MP2 calculations on a range of molecules, including linear chains and molecules with two- and three-dimensional character. The implementation is shared memory parallel via OpenMP and yields roughly 50% parallel efficiency with 16 cores for a large job. Considerable efforts were made to minimize memory demands, which increased as thresholds were tightened. A variety of relative energy calculations are presented as a function of threshold to provide some guidance to users on how to obtain adequate precision at a low compute cost. It is particularly clear that derivative properties require tighter thresholds in order to achieve an adequate precision.

Low-Scaling Algorithm for the Random Phase Approximation Using Tensor Hypercontraction with k-point Sampling

We present a low-scaling algorithm for the random phase approximation (RPA) with k-point sampling in the framework of tensor hypercontraction (THC) for electron repulsion integrals (ERIs). The THC factorization is obtained via a revised interpolative separable density fitting (ISDF) procedure with a momentum-dependent auxiliary basis for generic single-particle Bloch orbitals. Our formulation does not require preoptimized interpolating points or auxiliary bases, and the accuracy is systematically controlled by the number of interpolating points. The resulting RPA algorithm scales linearly with the number of k-points and cubically with the system size without any assumption on sparsity or locality of orbitals. The errors of ERIs and RPA energy show rapid convergence with respect to the size of the THC auxiliary basis, suggesting a promising and robust direction to construct efficient algorithms of higher order many-body perturbation theories for large-scale systems.

When Gold Is Not Enough: Platinum Standard of Quantum Chemistry with N7 Cost

In this paper, we extend the rank-reduced coupled-cluster formalism to the calculation of non-iterative energy corrections due to quadruple excitations. There are two major components of the proposed formalism. The first is an approximate compression of the quadruple excitation amplitudes using the Tucker format. The second is a modified functional used for the evaluation of the corrections which gives exactly the same results for the exact amplitudes, but is less susceptible to errors resulting from the aforementioned compression. We show, both theoretically and numerically, that the computational cost of the proposed method scales as the seventh power of the system size. Using reference results for a set of small molecules, the method is calibrated to deliver relative accuracy of a few percent in energy corrections. To illustrate the potential of the theory, we calculate the isomerization energy of ortho/meta benzyne (C₆H₄) and the barrier height for the Cope rearrangement in bullvalene (C₁₀H₁₀). The method retains a near-black-box nature of the conventional coupled-cluster formalism and depends on only one additional parameter that controls the accuracy.

Tensor-Hypercontracted MP2 First Derivatives: Runtime and Memory Efficient Computation of Hyperfine Coupling Constants

We employ our recently introduced tensor-hypercontracted (THC) second-order Møller-Plesset perturbation theory (MP2) method [Bangerter, F. H., Glasbrenner, M., Ochsenfeld, C. J. Chem. Theory Comput. 2021, 17, 211-221] for the computation of hyperfine coupling constants (HFCCs). The implementation leverages the tensor structure of the THC factorized electron repulsion integrals for an efficient formation of the integral-based intermediates. The computational complexity of the most expensive and formally quintic scaling exchange-like contribution is reduced to effectively subquadratic, by making use of the intrinsic, exponentially decaying coupling between tensor indices through screening based on natural blocking. Overall, this yields an effective subquadratic scaling with a low prefactor for the presented THC-based AO-MP2 method for the computation of isotropic HFCCs on DNA fragments with up to 500 atoms and 5000 basis functions. Furthermore, the implementation achieves considerable speedups with up to a factor of roughly 600-1000 compared to previous implementations [Vogler, S., Ludwig, M., Maurer, M., Ochsenfeld, C. J. Chem. Phys. 2017, 147, 024101] for medium-sized organic radicals, while also significantly reducing storage requirements.

Near-Exact CCSDT Energetics from Rank-Reduced Formalism Supplemented by Non-iterative Corrections

We introduce a non-iterative energy correction, added on top of the rank-reduced coupled-cluster method with single, double, and triple substitutions, that accounts for excitations excluded from the parent triple excitation subspace. The formula for the correction is derived by employing the coupled-cluster Lagrangian formalism, with an additional assumption that the parent excitation subspace is closed under the action of the Fock operator. Owing to the rank-reduced form of the triple excitation amplitudes tensor, the computational cost of evaluating the correction scales as N⁷, where N is the system size. The accuracy and computational efficiency of the proposed method is assessed for both total and relative correlation energies. We show that the non-iterative correction can fulfill two separate roles. If the accuracy level of a fraction of kJ/mol is sufficient for a given system, the correction significantly reduces the dimension of the parent triple excitation subspace needed in the iterative part of the calculations. Simultaneously, it enables reproducing the exact CCSDT results to an accuracy level below 0.1 kJ/mol, with a larger, yet still reasonable, dimension of the parent excitation subspace. This typically can be achieved at a computational cost only several times larger than required for the CCSD(T) method. The proposed method retains the black-box features of the single-reference coupled-cluster theory; the dimension of the parent excitation subspace remains the only additional parameter that has to be specified.

Hitting the Trifecta: How to Simultaneously Push the Limits of Schrödinger Solution with Respect to System Size, Convergence Accuracy, and Number of Computed States

Methods for solving the Schrödinger equation without approximation are in high demand but are notoriously computationally expensive. In practical terms, there are just three primary factors that currently limit what can be achieved: 1) system size/dimensionality; 2) energy level excitation; and 3) numerical convergence accuracy. Broadly speaking, current methods can deliver on any two of these three goals, but achieving all three at once remains an enormous challenge. In this paper, we shall demonstrate how to "hit the trifecta" in the context of molecular vibrational spectroscopy calculations. In particular, we compute the lowest 1000 vibrational states for the six-atom acetonitrile molecule (CH₃CN), to a numerical convergence of accuracy 10^-2 cm^-1 or better. These calculations encompass all vibrational states throughout most of the dynamically relevant range (i.e., up to ∼4250 cm^-1 above the ground state), computed in full quantum dimensionality (12 dimensions), to near spectroscopic accuracy. To our knowledge, no such vibrational spectroscopy calculation has ever previously been performed.

GPU acceleration of rank-reduced coupled-cluster singles and doubles

We have developed a graphical processing unit (GPU) accelerated implementation of our recently introduced rank-reduced coupled-cluster singles and doubles (RR-CCSD) method. RR-CCSD introduces a low-rank approximation of the doubles amplitudes. This is combined with a low-rank approximation of the electron repulsion integrals via Cholesky decomposition. The result of these two low-rank approximations is the replacement of the usual fourth-order CCSD tensors with products of second- and third-order tensors. In our implementation, only a single fourth-order tensor must be constructed as an intermediate during the solution of the amplitude equations. Owing in large part to the compression of the doubles amplitudes, the GPU-accelerated implementation shows excellent parallel efficiency (95% on eight GPUs). Our implementation can solve the RR-CCSD equations for up to 400 electrons and 1550 basis functions-roughly 50% larger than the largest canonical CCSD computations that have been performed on any hardware. In addition to increased scalability, the RR-CCSD computations are faster than the corresponding CCSD computations for all but the smallest molecules. We test the accuracy of RR-CCSD for a variety of chemical systems including up to 1000 basis functions and determine that accuracy to better than 0.1% error in the correlation energy can be achieved with roughly 95% compression of the ov space for the largest systems considered. We also demonstrate that conformational energies can be predicted to be within 0.1 kcal mol^-1 with efficient compression applied to the wavefunction. Finally, we find that low-rank approximations of the CCSD doubles amplitudes used in the similarity transformation of the Hamiltonian prior to a conventional equation-of-motion CCSD computation will not introduce significant errors (on the order of a few hundredths of an electronvolt) into the resulting excitation energies.

A diagrammatic approach for automatically deriving analytical gradients of tensor hyper-contracted electronic structure methods

We introduce a diagrammatic approach to facilitate the automatic derivation of analytical nuclear gradients for tensor hyper-contraction (THC) based electronic structure methods. The automatically derived gradients are guaranteed to have the same scaling in terms of both operation count and memory footprint as the underlying energy calculations, and the computation of a gradient is roughly three times as costly as the underlying energy. The new diagrammatic approach enables the first cubic scaling implementation of nuclear derivatives for THC tensors fitted in molecular orbital basis (MO-THC). Furthermore, application of this new approach to THC-MP2 analytical gradients leads to an implementation, which is at least four times faster than the previously reported, manually derived implementation. Finally, we apply the new approach to the 14 tensor contraction patterns appearing in the supporting subspace formulation of multireference perturbation theory, laying the foundation for developments of analytical nuclear gradients and nonadiabatic coupling vectors for multi-state CASPT2.

Cubic-Scaling All-Electron GW Calculations with a Separable Density-Fitting Space-Time Approach

We present an implementation of the GW space-time approach that allows cubic-scaling all-electron calculations with standard Gaussian basis sets without exploiting any localization or sparsity considerations. The independent-electron susceptibility is constructed in a time representation over a nonuniform distribution of real-space locations {r_k} optimized within a separable resolution-of-the-identity framework to reproduce standard Coulomb-fitting calculations with meV accuracy. The compactness of the obtained {r_k} distribution leads to a crossover with the standard Coulomb-fitting scheme for system sizes below a few hundred electrons. The needed analytic continuation follows a recent approach that requires the continuation of the screened Coulomb potential rather than the much more structured self-energy. The present scheme is benchmarked over large molecular sets, and scaling properties are demonstrated on a family of defected hexagonal boron-nitride flakes containing up to 6000 electrons.

A critical analysis of least-squares tensor hypercontraction applied to MP3

The least-squares tensor hypercontraction (LS-THC) approach is a promising method of reducing the high polynomial scaling of wavefunction methods, for example, those based on many-body perturbation theory or coupled cluster. Here, we focus on LS-THC-MP3 and identify four variants with differing errors and efficiency characteristics. The performance of LS-THC-MP3 is analyzed for regular test systems with up to 40 first-row atoms. We also analyze the size-extensivity/size-consistency and grid- and basis set dependence of LS-THC-MP3. Overall, the errors observed are favorably small in comparison with standard density fitting, and a more streamlined method of generating grids via pruning is suggested. A practical crossover (the point at which LS-THC-MP3 is cheaper than the canonical method) is achieved around 240 correlated electrons. Despite several drawbacks of LS-THC that have been identified: an initial non-linearity of error when increasing system size, poor description of angular correlation, and a potentially large increase in error with the basis set size, the results show that LS-THC has significant potential for practical application to MP3 and other wavefunction methods.

Reduced scaling formulation of CASPT2 analytical gradients using the supporting subspace method

We present a reduced scaling and exact reformulation of state specific complete active space second-order perturbation (CASPT2) analytical gradients in terms of the MP2 and Fock derivatives using the supporting subspace method. This work follows naturally from the supporting subspace formulation of the CASPT2 energy in terms of the MP2 energy using dressed orbitals and Fock builds. For a given active space configuration, the terms corresponding to the MP2-gradient can be evaluated with O(N⁵) operations, while the rest of the calculations can be computed with O(N³) operations using Fock builds, Fock gradients, and linear algebra. When tensor-hyper-contraction is applied simultaneously, the computational cost can be further reduced to O(N⁴) for a fixed active space size. The new formulation enables efficient implementation of CASPT2 analytical gradients by leveraging the existing graphical processing unit (GPU)-based MP2 and Fock routines. We present benchmark results that demonstrate the accuracy and performance of the new method. Example applications of the new method in ab initio molecular dynamics simulation and constrained geometry optimization are given.

Low-Scaling Tensor Hypercontraction in the Cholesky Molecular Orbital Basis Applied to Second-Order Møller-Plesset Perturbation Theory

We employ various reduced scaling techniques to accelerate the recently developed least-squares tensor hypercontraction (LS-THC) approximation [Parrish, R. M., Hohenstein, E. G., Martínez, T. J., Sherrill, C. D. J. Chem. Phys. 137, 224106 (2012)] for electron repulsion integrals (ERIs) and apply it to second-order Møller-Plesset perturbation theory (MP2). The grid-projected ERI tensors are efficiently constructed using a localized Cholesky molecular orbital basis from density-fitted integrals with an attenuated Coulomb metric. Additionally, rigorous integral screening and the natural blocking matrix format are applied to reduce the complexity of this step. By recasting the equations to form the quantized representation of the 1/r operator Z into the form of a system of linear equations, the bottleneck of inverting the grid metric via pseudoinversion is removed. This leads to a reduced scaling THC algorithm and application to MP2 yields the (sub-)quadratically scaling THC-ω-RI-CDD-SOS-MP2 method. The efficiency of this method is assessed for various systems including DNA fragments with over 8000 basis functions and the subquadratic scaling is illustrated.

Reduced scaling extended multi-state CASPT2 (XMS-CASPT2) using supporting subspaces and tensor hyper-contraction

We present a reduced scaling formulation of the extended multi-state CASPT2 (XMS-CASPT2) method, which is based on our recently developed state-specific CASPT2 (SS-CASPT2) formulation using supporting subspaces and tensor hyper-contraction. By using these two techniques, the off-diagonal elements of the effective Hamiltonian can be computed with only O(N³) operations and O(N²) memory, where N is the number of basis functions. This limits the overall computational scaling to O(N⁴) operations and O(N²) memory. Thus, excited states can now be obtained at the same reduced (relative to previous algorithms) scaling we achieved for SS-CASPT2. In addition, we also investigate how the energy denominators can be factorized with the Laplace quadrature when some of the denominators are negative, which is critical for excited state calculations. An efficient implementation of the method has been developed using graphical processing units while also exploiting spatial sparsity in tensor operations. We benchmark the accuracy of the new method by comparison to non-THC formulated XMS-CASPT2 for the excited states of various molecules. In our tests, the THC approximation introduces negligible errors (≈0.01 eV) compared to the non-THC reference method. Scaling behavior and computational timings are presented to demonstrate performance. The new method is also interfaced with quantum mechanics/molecular mechanics (QM/MM). In an example study of green fluorescent protein, we show how the XMS-CASPT2 potential energy surfaces and excitation energies are affected by increasing the size of the QM region up to 278 QM atoms with more than 2300 basis functions.

TeraChem: Accelerating electronic structure and ab initio molecular dynamics with graphical processing units

Developed over the past decade, TeraChem is an electronic structure and ab initio molecular dynamics software package designed from the ground up to leverage graphics processing units (GPUs) to perform large-scale ground and excited state quantum chemistry calculations in the gas and the condensed phase. TeraChem's speed stems from the reformulation of conventional electronic structure theories in terms of a set of individually optimized high-performance electronic structure operations (e.g., Coulomb and exchange matrix builds, one- and two-particle density matrix builds) and rank-reduction techniques (e.g., tensor hypercontraction). Recent efforts have encapsulated these core operations and provided language-agnostic interfaces. This greatly increases the accessibility and flexibility of TeraChem as a platform to develop new electronic structure methods on GPUs and provides clear optimization targets for emerging parallel computing architectures.

A Quadratic Pair Atomic Resolution of the Identity Based SOS-AO-MP2 Algorithm Using Slater Type Orbitals

We report a production level implementation of pair atomic resolution of the identity (PARI) based second-order Møller-Plesset perturbation theory (MP2) in the Slater type orbital (STO) based Amsterdam Density Functional (ADF) code. As demonstrated by systematic benchmarks, dimerization and isomerization energies obtained with our code using STO basis sets of triple-ζ-quality show mean absolute deviations from Gaussian type orbital, canonical, basis set limit extrapolated, global density fitting (DF)-MP2 results of less than 1 kcal/mol. Furthermore, we introduce a quadratic scaling atomic orbital based spin-opposite-scaled (SOS)-MP2 approach with a very small prefactor. Due to a worst-case scaling of [Formula: see text], our implementation is very fast already for small systems and shows an exceptionally early crossover to canonical SOS-PARI-MP2. We report computational wall time results for linear as well as for realistic three-dimensional molecules and show that triple-ζ quality calculations on molecules of several hundreds of atoms are only a matter of a few hours on a single compute node, the bottleneck of the computations being the SCF rather than the post-SCF energy correction.

Q-MP2-OS: Møller-Plesset Correlation Energy by Quadrature

We present a quadrature-based algorithm for computing the opposite-spin component of the MP2 correlation energy which scales quadratically with basis set size and is well-suited to large-scale parallelization. The key ideas, which are rooted in the earlier work of Hirata and co-workers, are to abandon all two-electron integrals, recast the energy as a seven-dimensional integral, approximate that integral by quadrature, and employ a cutoff strategy to minimize the number of intermediate quantities. We discuss our implementation in detail and show that it parallelizes almost perfectly on 840 cores for cyclosporine (a molecule with roughly 200 atoms), exhibits [Formula: see text] scaling for a sequence of polyglycines, and is principally limited by the accuracy of its quadrature.

A direct-compatible formulation of the coupled perturbed complete active space self-consistent field equations on graphical processing units

We recently developed an algorithm to compute response properties for the state-averaged complete active space self-consistent field method (SA-CASSCF) that capitalized on sparsity in the atomic orbital basis. Our original algorithm was limited to treating small to moderate sized active spaces, but the recent development of graphical processing unit (GPU) based direct-configuration interaction algorithms provides an opportunity to extend this to large active spaces. We present here a direct-compatible version of the coupled perturbed equations, enabling us to compute response properties for systems treated with arbitrary active spaces (subject to available memory and computation time). This work demonstrates that the computationally demanding portions of the SA-CASSCF method can be formulated in terms of seven fundamental operations, including Coulomb and exchange matrix builds and their derivatives, as well as, generalized one- and two-particle density matrix and σ vector constructions. As in our previous work, this algorithm exhibits low computational scaling and is accelerated by the use of GPUs, making possible optimizations and nonadiabatic dynamics on systems with O(1000) basis functions and O(100) atoms, respectively.

The Laplace transformed divide-expand-consolidate resolution of the identity second-order Møller-Plesset perturbation (DEC-LT-RIMP2) theory method

The divide-expand-consolidate resolution of the identity second-order Møller-Plesset perturbation (DEC-RI-MP2) theory method introduced in Baudin et al. [J. Chem. Phys. 144, 054102 (2016)] is significantly improved by introducing the Laplace transform of the orbital energy denominator in order to construct the double amplitudes directly in the local basis. Furthermore, this paper introduces the auxiliary reduction procedure, which reduces the set of the auxiliary functions employed in the individual fragments. The resulting Laplace transformed divide-expand-consolidate resolution of the identity second-order Møller-Plesset perturbation method is applied to the insulin molecule where we obtain a factor 9.5 speedup compared to the DEC-RI-MP2 method.