The divide–expand–consolidate MP2 scheme goes massively parallel

An Efficient RI-MP2 Algorithm for Distributed Many-GPU Architectures

Second-order Møller-Plesset perturbation theory (MP2) using the Resolution of the Identity approximation (RI-MP2) is a widely used method for computing molecular energies beyond the Hartree-Fock mean-field approximation. However, its high computational cost and lack of efficient algorithms for modern supercomputing architectures limit its applicability to large molecules. In this paper, we present the first distributed-memory many-GPU RI-MP2 algorithm explicitly designed to utilize hundreds of GPU accelerators for every step of the computation. Our novel algorithm achieves near-peak performance on GPU-based supercomputers through the development of a distributed memory algorithm for forming RI-MP2 intermediate tensors with zero internode communication, except for a single O ( N 2 ) asynchronous broadcast, and a distributed memory algorithm for the O ( N 5 ) energy reduction step, capable of sustaining near-peak performance on clusters with several hundred GPUs. Comparative analysis shows our implementation outperforms state-of-the-art quantum chemistry software by over 3.5 times in speed while achieving an 8-fold reduction in computational power consumption. Benchmarking on the Perlmutter supercomputer, our algorithm achieves 11.8 PFLOP/s (83% of peak performance) performing and the RI-MP2 energy calculation on a 314-water cluster with 7850 primary and 30,144 auxiliary basis functions in 4 min on 180 nodes and 720 A100 GPUs. This performance represents a substantial improvement over traditional CPU-based methods, demonstrating significant time-to-solution and power consumption benefits of leveraging modern GPU-accelerated computing environments for quantum chemistry calculations.

Local Second Order Mo̷ller-Plesset Theory with a Single Threshold Using Orthogonal Virtual Orbitals: A Distributed Memory Implementation

In order to alleviate the computational burden associated with superlinear compute scalings with molecular size in electron correlation methods, researchers have developed local correlation methods that wisely treat relatively small contributions as zeros but still yield accurate energy approximation. Such local correlation techniques can also be combined with parallel computing resources to obtain further efficiency and scalability. This work focuses on the distributed memory parallel implementation of a local correlation method for second order Mo̷ller-Plesset (MP2) theory. This method also only has a single threshold to control the dropping of terms and accuracy of different computing kernels in the algorithm. The process partitioning strategy and distributed parallel implementation with the message passing interface (MPI) are discussed. In particular, the algorithm relies on a fixed sparsity pattern matrix multiplication and a corresponding distributed conjugate gradient solver, which exhibits almost linear scaling in both strong and weak scaling analyses. Numerical experiments on a range of molecules, including linear chains and molecules with 2 and 3-dimensional characters, are reported. For example, with only 32 MPI ranks, this MP2 implementation can calculate the correlation energy of vancomycin in def2-TZVP basis within 0.003% accuracy (10-6.5 threshold) in half an hour, where the same problem is unfeasible to solve with sequential or pure shared memory implementations.

Accurate Calculation of Isomerization and Conformational Energies of Larger Molecules Using Explicitly Correlated Local Coupled Cluster Methods in Molpro and ORCA

An overview of the approximations in the explicitly correlated local coupled cluster methods PNO-LCCSD(T)-F12 in Molpro and DLPNO-CCSD(T)F12 in ORCA is given. Options to select the domains of projected atomic orbitals (PAOs), pair natural orbitals (PNOs), and triples natural orbitals (TNOs) in both programs are described and compared in detail. The two programs are applied to compute isomerization and conformational energies of the ISOL24 and ACONFL test sets, where the former is part of the GMTKN55 benchmark suite. Thorough studies of basis set effects are presented for selected systems. These revealed large intramolecular basis set superposition effects that make it practically impossible to reliably determine the complete basis set (CBS) limits without including explicitly correlated terms. The latter strongly reduce the basis set dependence and at the same time also errors caused by the local domain approximations. On the basis of these studies, the PNO-LCCSD(T)-F12 method is applied to determine new reference energies for the above-mentioned benchmark sets. We are confident that our results should agree within a few tenths of a kcal mol-1 with the (unknown) CCSD(T)/CBS values, which therefore allowed us to define computational settings for accurate explicitly correlated local coupled cluster methods with moderate computational effort. With these protocols, especially PNO-LCCSD(T)-F12b/AVTZ', reliable reference values for comprehensive benchmark sets can be generated efficiently. This can significantly advance the development and evaluation of the performance of approximate electronic structure methods, especially improved density functional approximations or machine learning approaches.

Direct orbital selection within the domain-based local pair natural orbital coupled-cluster method

Domain-based local pair natural orbital coupled cluster (DLPNO-CC) has become increasingly popular to calculate relative energies (e.g., reaction energies and reaction barriers). It can be applied within a multi-level DLPNO-CC-in-DLPNO-CC ansatz to reduce the computational cost and focus the available computational resources on a specific subset of the occupied orbitals. We demonstrate how this multi-level DLPNO-CC ansatz can be combined with our direct orbital selection (DOS) approach [M. Bensberg and J. Neugebauer, J. Chem. Phys. 150, 214106 (2019)] to automatically select orbital sets for any multi-level calculation. We find that the parameters for the DOS procedure can be chosen conservatively such that they are transferable between reactions. The resulting automatic multi-level DLPNO-CC method requires no user input and is extremely robust and accurate. The computational cost is easily reduced by a factor of 3 without sacrificing accuracy. We demonstrate the accuracy of the method for a total of 61 reactions containing up to 174 atoms and use it to predict the relative stability of conformers of a Ru-based catalyst.

Scalable Electron Correlation Methods. 8. Explicitly Correlated Open-Shell Coupled-Cluster with Pair Natural Orbitals PNO-RCCSD(T)-F12 and PNO-UCCSD(T)-F12

We present explicitly correlated open-shell pair natural orbital local coupled-cluster methods, PNO-RCCSD(T)-F12 and PNO-UCCSD(T)-F12. The methods are extensions of our previously reported PNO-R/UCCSD methods (J. Chem. Theory Comput., 2020, 16, 3135-3151, https://pubs.acs.org/doi/10.1021/acs.jctc.0c00192) with additions of explicit correlation and perturbative triples corrections. The explicit correlation treatment follows the spin-orbital CCSD-F12b theory using Ansatz 3*A, which is found to yield comparable or better basis set convergence than the more rigorous Ansatz 3C in computed ionization potentials and reaction energies using double- to quaduple-ζ basis sets. The perturbative triples correction is adapted from the spin-orbital (T) theory to use triples natural orbitals (TNOs). To address the coupling due to off-diagonal Fock matrix elements, the local triples amplitudes are iteratively solved using small domains of TNOs, and a semicanonical (T0) domain correction with larger domains is applied to reduce the domain errors. The performance of the methods is demonstrated through benchmark calculations on ionization potentials, radical stabilization energies, reaction energies of fragmentations and rearrangements in radical cations, and spin-state energy differences of iron complexes. For a few test sets where canonical calculations are feasible, PNO-RCCSD(T)-F12 results agree with the canonical ones to within 0.4 kcal mol^-1, and this maximum error is reduced to below 0.2 kcal mol^-1 when large local domains are used. For larger systems, results using different thresholds for the local approximations are compared to demonstrate that 1 kcal mol^-1 level of accuracy can be achieved using our default settings. For a couple of difficult cases, it is demonstrated that the errors from individual approximations are only a fraction of 1 kcal mol^-1, and the overall accuracy of the method does not rely on error compensations. In contrast to canonical calculations, the use of spin-orbitals does not lead to a significant increase of computational time and memory usage in the most expensive steps of PNO-R/UCCSD(T)-F12 calculations. The only exception is the iterative solution of the (T) amplitudes, which can be avoided without significant errors by using a perturbative treatment of the off-diagonal coupling, known as (T1) approximation. For most systems, even the semicanonical approximation (T0) leads only to small errors in relative energies. Our program is well parallelized and capable of computing accurate correlation energies for molecules with 100-200 atoms using augmented triple-ζ basis sets in less than a day of elapsed time on a small computer cluster.

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.

Implementation of the iterative triples model CC3 for excitation energies using pair natural orbitals and Laplace transformation techniques

We present a pair natural orbital (PNO)-based implementation of CC3 excitation energies, which extends our previously published state-specific PNO ansatz for the solution of the excited state eigenvalue problem to methods including connected triple excitations. A thorough analysis of the equations for the excited state triples amplitudes is presented from which we derive a suitable state-specific triple natural orbital basis for the excited state triples amplitudes, which performs equally well for local and non-local excitations. The accuracy of the implementation is evaluated using a large and diverse test set. We find that for states with small contributions from double excitations, a T0 approximation to PNO-CC3 yields accurate results with a mean absolute error (MAE) for T_PNO = 10^-7 in the range of 0.02 eV. However, for states with larger double excitation contributions, the T0 approximation is found to yield significantly less accurate results, while the Laplace-transformed variant of PNO-CC3 shows a uniform accuracy for singly and doubly excited states (MAE and maximum error of 0.01 eV and 0.07 eV for T_PNO = 10^-7, respectively). Finally, we apply PNO-CC3 to the calculation of the first excited state of berenil at a S₁ minimum geometry, which is shown to be close to a conical intersection. This calculation in the aug-cc-pVTZ basis set (more than 1300 basis functions) is the largest calculation ever performed with CC3 on excitation energies.

Double hybrid DFT calculations with Slater type orbitals

On a comprehensive database with 1,644 datapoints, covering several aspects of main-group as well as of transition metal chemistry, we assess the performance of 60 density functional approximations (DFA), among them 36 double hybrids (DH). All calculations are performed using a Slater type orbital (STO) basis set of triple-ζ (TZ) quality and the highly efficient pair atomic resolution of the identity approach for the exchange- and Coulomb-term of the KS matrix (PARI-K and PARI-J, respectively) and for the evaluation of the MP2 energy correction (PARI-MP2). Employing the quadratic scaling SOS-AO-PARI-MP2 algorithm, DHs based on the spin-opposite-scaled (SOS) MP2 approximation are benchmarked against a database of large molecules. We evaluate the accuracy of STO/PARI calculations for B3LYP as well as for the DH B2GP-PLYP and show that the combined basis set and PARI-error is comparable to the one obtained using the well-known def2-TZVPP Gaussian-type basis set in conjunction with global density fitting. While quadruple-ζ (QZ) calculations are currently not feasible for PARI-MP2 due to numerical issues, we show that, on the TZ level, Jacob's ladder for classifying DFAs is reproduced. However, while the best DHs are more accurate than the best hybrids, the improvements are less pronounced than the ones commonly found on the QZ level. For conformers of organic molecules and noncovalent interactions where very high accuracy is required for qualitatively correct results, DHs provide only small improvements over hybrids, while they still excel in thermochemistry, kinetics, transition metal chemistry and the description of strained organic systems.

The Quest for Highly Accurate Excitation Energies: A Computational Perspective

We provide an overview of the successive steps that made it possible to obtain increasingly accurate excitation energies with computational chemistry tools, eventually leading to chemically accurate vertical transition energies for small- and medium-size molecules. First, we describe the evolution of ab initio methods employed to define benchmark values, with the original Roos CASPT2 method, then the CC3 method as in the renowned Thiel set, and more recently the resurgence of selected configuration interaction methods. The latter method has been able to deliver consistently, for both single and double excitations, highly accurate excitation energies for small molecules, as well as medium-size molecules with compact basis sets. Second, we describe how these high-level methods and the creation of representative benchmark sets of excitation energies have allowed the fair and accurate assessment of the performance of computationally lighter methods. We conclude by discussing possible future theoretical and technological developments in the field.

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.

Scalable Electron Correlation Methods. 6. Local Spin-Restricted Open-Shell Second-Order Møller-Plesset Perturbation Theory Using Pair Natural Orbitals: PNO-RMP2

We present a (near) linear scaling implementation of high-spin open-shell Møller-Plesset perturbation theory using pair natural orbitals (PNO-RMP2). The theory is based on a new variant of open-shell MP2 which is fully spin-adapted and uses a single set of spin-free amplitudes, as in closed-shell MP2. This method, denoted SROMP2, is invariant to unitary orbital transformations within the closed, open, and virtual orbital subspaces. Accordingly, only a single set of PNOs per spatial orbital pair is needed, and the efficiency is similar to closed-shell calculations. The PNOs are obtained using a semicanonical approximation with large domains of projected atomic orbitals (PAOs). Linear scaling is achieved provided that the open-shell orbitals are local, and distant pairs are treated by multipole approximations. The method is efficiently parallelized. The convergence of ionization and reaction energies as a function of the PAO and PNO domain sizes is demonstrated and found to be very similar as for closed-shell calculations. The suitability of the PNOs for explicitly correlated PNO-RCCSD-F12 calculations is also tested. So far, this method is only simulated using a conventional program with appropriate projections to the PAO and PNO subspaces. It is demonstrated for radical stabilization energies as well as ionization potentials that the errors caused by the local domain approximations with our default thresholds are negligible.

Low-scaling first-order properties within second-order Møller-Plesset perturbation theory using Cholesky decomposed density matrices

An efficient implementation of energy gradients and of hyperfine coupling constants in second-order Møller-Plesset perturbation theory (MP2) is presented based on our fully atomic orbital (AO)-based approach. For the latter, an unrestricted AO-based MP2 formulation is introduced. A reduction in the dependency of the computational efficiency on the size of the basis set is achieved by a Cholesky decomposition and the prefactor is reduced by the resolution-of-the-identity approximation. Significant integral contributions are selected based on distance-including integral estimates (denoted as QQR-screening) and its reliability as a fully controlled screening procedure is demonstrated. The rate-determining steps are shown via model computations to scale cubically in the computation of energy gradients and quadratically in the case of hyperfine coupling constants. Furthermore, a significant speed-up of the computational time with respect to the canonical formulation is demonstrated.

Divide-Expand-Consolidate Second-Order Møller-Plesset Theory with Periodic Boundary Conditions

We present a generalization of the divide-expand-consolidate (DEC) framework for local coupled-cluster calculations to periodic systems and test it at the second-order Møller-Plesset (MP2) level of theory. For simple model systems with periodicity in one, two, and three dimensions, comparisons with extrapolated molecular calculations and the local MP2 implementation in the Cryscor program show that the correlation energy errors of the extended DEC (X-DEC) algorithm can be controlled through a single parameter, the fragment optimization threshold. Two computational bottlenecks are identified: the size of the virtual orbital spaces and the number of pair fragments required to achieve a given accuracy of the correlation energy. For the latter, we propose an affordable algorithm based on cubic splines interpolation of a limited number of pair-fragment interaction energies to determine a pair cutoff distance in accordance with the specified fragment optimization threshold.

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.

Efficient linear-scaling second-order Møller-Plesset perturbation theory: The divide-expand-consolidate RI-MP2 model

The Resolution of the Identity second-order Møller-Plesset perturbation theory (RI-MP2) method is implemented within the linear-scaling Divide-Expand-Consolidate (DEC) framework. In a DEC calculation, the full molecular correlated calculation is replaced by a set of independent fragment calculations each using a subset of the total orbital space. The number of independent fragment calculations scales linearly with the system size, rendering the method linear-scaling and massively parallel. The DEC-RI-MP2 method can be viewed as an approximation to the DEC-MP2 method where the RI approximation is utilized in each fragment calculation. The individual fragment calculations scale with the fifth power of the fragment size for both methods. However, the DEC-RI-MP2 method has a reduced prefactor compared to DEC-MP2 and is well-suited for implementation on massively parallel supercomputers, as demonstrated by test calculations on a set of medium-sized molecules. The DEC error control ensures that the standard RI-MP2 energy can be obtained to the predefined precision. The errors associated with the RI and DEC approximations are compared, and it is shown that the DEC-RI-MP2 method can be applied to systems far beyond the ones that can be treated with a conventional RI-MP2 implementation.

First UHF Implementation of the Incremental Scheme for Open-Shell Systems

The incremental scheme makes it possible to compute CCSD(T) correlation energies to high accuracy for large systems. We present the first extension of this fully automated black-box approach to open-shell systems using an Unrestricted Hartree-Fock (UHF) wave function, extending the efficient domain-specific basis set approach to handle open-shell references. We test our approach on a set of organic and metal organic structures and molecular clusters and demonstrate standard deviations from canonical CCSD(T) values of only 1.35 kJ/mol using a triple ζ basis set. We find that the incremental scheme is significantly more cost-effective than the canonical implementation even for relatively small systems and that the ease of parallelization makes it possible to perform high-level calculations on large systems in a few hours on inexpensive computers. We show that the approximations that make our approach widely applicable are significantly smaller than both the basis set incompleteness error and the intrinsic error of the CCSD(T) method, and we further demonstrate that incremental energies can be reliably used in extrapolation schemes to obtain near complete basis set limit CCSD(T) reaction energies for large systems.