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.

Introduction

Spurred by the interest in large biomolecules and inorganic systems, the last decades have witnessed a tremendous effort in making accurate electronic structure methods routinely applicable to molecules and solids of ever increasing size. Due to its still unrivaled price/performance-ratio,[1] Kohn–Sham (KS)[2] density functional theory (DFT)[3,4] has established itself as the workhorse of quantum chemistry for medium and large systems.[5−11] Unfortunately, due to an insufficient description of electron correlation, state of the art semilocal[12] or hybrid[13] approximations to the exact exchange–correlation functional often fail to accurately account for London dispersion-type effects[14−18] and noncovalent interactions.[19] Both are of paramount importance for a thorough understanding of the properties and reactivity of biochemical systems and organometallic compounds.[20,21]

Wave function based ab initio methods, however, offer a systematic route toward the complete and explicit description of all dynamical correlation effects. As known only too well, this does not come for free: Expressed in terms of canonical orbitals, their steep computational scaling (N5 (N is a measure for the system size[22]) for second-order Møller–Plesset perturbation theory (MP2),[23]N6 for CC with singles and doubles excitations (CCSD),[24−26]N7 for CCSD with a perturbative treatment of triple excitations (CCSD(T)),[27] respectively), their tremendous memory requirements in practical implementations,[28] and their slow convergence to the basis set-limit[29] complicate the application of these methods to large molecular systems.

The commonly realized concept of the locality of dynamical electron correlation has led to a family of low-scaling wave function based methods approaching the accuracy of their canonical counterparts. The field came to life when Pulay, Saebø, and co-workers employed localized molecular orbitals (MO)[30−36] and restricted excitations to local domains in configuration interaction (CI)[37,38] and MPn[39−41] computations. Successfully transferred to the realms of highly accurate CC theory by Werner and co-workers,[42] a plethora of low-scaling CC[43−49] and MP2[51−55] codes has been developed. With the size of the excitation domains becoming a liming factor, Neese and co-workers[28,56−59] and others[60−69] revived the decades ago developed[70−83] pair natural orbital (PNO)[84,85] approach to further compress the virtual subspaces. At the same time, relying on a physical partition of the system of interest instead of the orbital space, several flavors of fragment based approaches have been brought forward. Most notably these are the incremental method,[86−88] the cluster-in-molecule,[89−91] the divide-and-conquer,[92−98] and the divide-expand-consolidate[99−106] approaches.

While having greatly extended the range of computationally tractable molecules,[58,68,69,107] local CC methods can still not compete with the wall clock times of KS-DFT computations on the hybrid level. This is not necessarily the case for MP2.[59,65,67] Although less accurate than CCSD and higher CC methods, MP2 has been demonstrated to accurately describe properties such as dispersion interactions and hydrogen binding,[18,108,109] NMR chemical shifts,[110] and polarizabilities[111,112] as well as molecular interaction energies[113] and conformational energies,[114] especially in its spin-scaled variants.[115−118] Furthermore, MP2 energies from KS orbitals have been used extensively to incorporate explicit electron correlation effects into the calculation of KS-DFT energies.[119−125] These so-called double hybrid (DH) density functionals often yield clear improvements over hybrid functionals,[126,127] making low-cost MP2 implementations extremely desirable.

While a simple expression for the MP2 amplitudes exists using canonical orbitals, the MP2 amplitude equations can only be solved iteratively in most fragments, as well as localized orbital based formalisms (the integral-direct formulation of Nagy et al.[67,91] should be mentioned as a prominent exception). Retaining the simplicity of the canonical formalism, exploiting sparsity in the electron-repulsion integral (ERI) tensor[128] instead is yet another alley toward the order-N computation of MP2 energies. Formally quartically scaling, it is well-known that the number of nonzero elements in the ERI tensor only scales linearly with system size. The resolution of the identity (RI) or density fitting approximation (DF)[129−140] is the most popular technique to overcome the scaling of the fourth-order ERI tensor by decomposing it into third-order and second-order tensors and has been applied successfully to reduce the prefactor of canonical MP2.[48,105,106,141−149]

Reformulating the energy denominator of canonical MP2 as an integral expression[150−152] (often referred to as a Laplace transform (LT)), the MP2 energy can be evaluated in the atomic orbital (AO) basis. In this representation, approaches to reduce the dimensionality of the ERI tensor can be applied more efficiently. Using integral prescreening techniques,[152−156] reduced-scaling AO-MP2 codes could be realized by Ayala and Scuseria,[157,158] Ochsenfeld and co-workers,[50,159−162] and others.[163] However, employing large sets of spatially extended AOs, the prefactor of these methods is increased significantly due to an increase of linear dependencies in the pseudo density matrices (PDM).[160,164] Cholesky decomposition (CD) based techniques[165−169] can be used to obtain a set of localized occupied (virtual) orbitals with cardinality equal to the rank of the occupied (virtual) PDMs. Employing this approach in conjunction with screening techniques and DF, very fast AO-MP2 implementations have been reported by the Ochsenfeld group.[164,170−173] Also the tensor hypercontraction (THC) approach to compress the ERI tensor, brought forward by Martínez and co-workers,[174−178] has been used successfully to obtain fast and reliable low-scaling MP2 implementations.[174,179−181]

In this work, we explore the use of the pair atomic resolution of the identity (PARI)[129,182−187] in the evaluation of MP2 energies. In this approach, the cubic scaling of global DF is reduced automatically by expanding products of AO pairs in terms of fit functions centered on the same atoms as the pair of target functions only. In this way, the complexity of the ERI tensor only scales quadratically with system size and linearly when insignificant pairs of basis functions are neglected. Already developed as early as in the 1970s by Baerends et al.,[129] its potential to accelerate the evaluation of the exact exchange has only been demonstrated recently.[183,185−188] Although some concerns regarding its accuracy and numerical stability have been brought forward,[183,189] the applicability of the PARI approximation in the framework of MP2 and the random phase approximation (RPA) has been reported by Ihrig et al.[186]

The purpose of this paper is 2-fold. First, we present a production level implementation of PARI based MP2 (PARI-MP2) using STOs. By means of systematic benchmarks on various test sets, we demonstrate that our PARI-MP2 implementation in conjunction with triple-ζ quality basis set[190] reproduces the basis set limit of canonical MP2 within less than 1 kcal/mol on average. Second, we show how the PARI approach can be used to considerably speed up the evaluation of MP2 energies in a memory efficient way. As we think that these are the most realistic targets for local correlation methods, we explicitly focused on designing an algorithm which is very fast and reliable for small as well as medium, compact systems of up to several hundreds of atoms: We herein report a quadratic scaling AO based spin-opposite-scaled (SOS)-PARI-MP2 implementation of order-N3 without any distance effects considered. Due to a small prefactor, the computation of the SOS-MP2 energy is much faster than the SCF itself. Although we herein refer often to STO basis sets, we emphasize that our approach is completely general and can be implemented with basis sets of arbitrary type.

This paper is organized as follows: In section , we introduce the PARI approach to evaluate the exact exchange and describe our AO based SOS-PARI-MP2 implementation (SOS-AO-PARI-MP2 in the following). In section , we report our benchmark results and demonstrate the excellent scaling properties of our algorithm for selected medium and large molecules. Finally, in section we summarize and conclude this work.

The herein frequently appearing indices are summarized as follows:

A, B, C, D: indices denoting atomic centers, where A = 1, ..., Na and Na denotes the number of atoms in the system.

μ, ν, κ, λ: AO basis indices, where μ = 1, ..., NbA, where NbA denotes the number of basis functions Nb centered on atom A.

α, β, γ, δ: auxiliary basis indices, where α = 1, ..., NfA, and NfA is the number of auxiliary fit functions Nf centered on atom A.

The following convention is always applied: μ, α ∈ A, ν, β ∈ B, κ, γ ∈ C, λ, δ ∈ D, which means that the indices μ,α always imply that the corresponding functions are centered on atom A.

μ̃ (α̃) denote global AO basis indices (global auxiliary basis indices), ranging from 1 to Nb,all (Nf,all), where Nb,all (N) are the number of AOs (auxiliary fit functions) of the whole system.

τ: Numerical quadrature point indices, ranging from 1 to Nnq, where Nnq denotes the number of quadrature points.

Theory

Pair Atomic Resolution of the Identity for the Exact Exchange

PARI (or local pair fitting) is a quite extreme variant of the more general local domain fitting approach. The general idea is to overcome the cubic scaling of global, Coulomb metric based DF using fit functions in the neighborhood of the target product only. When pair-atomic densities are chosen as local fitting domains, the PARI method is obtained. This approach is clearly more approximate than global DF but can be physically motivated as atom-centered basis functions generate strongly localized contributions to the total density; i.e., the AO pair density products χμ(r) χν(r) ≡ ρ are local in nature. Here, we only recapitulate shortly the concepts which are necessary to understand our MP2 implementation. We are planning to elaborate on our PARI-HF implementation in a forthcoming publication.

Each atomic pair density with AOs centered on a certain pair of atoms (A,B) is expanded in a local set of fit functions {fα(r), fβ(r)},where (in accordance with the convention introduced in section ) the fit functions are defined to be centered on the second atom for each pair product. Similarly to the global “RI-V” approach[132,136] we define a residual vectorand minimize the self-repulsion of the residual, (ε|ε). Essentially, this procedure minimizes the error in the electric field generated by the two charge distributions[191] and consequently minimizes the error in the representation of the ERIs.[136] It leads to a set of linear equations for the determination of the fit coefficients (where the inverse of the matrix V can always be calculated in a numerically stable way using a singular value decomposition (SVD)),with the electrostatic interaction between the product of AOs and fit functions,and the electrostatic interaction (Coulomb overlap) between fit functions,where the electrostatic potential v due to the function fα has been introduced as

Using (1) and (5), the elements of the ERI tensor (in Mulliken notation) might be expressed asWith the set of eqs –7, we have essentially expressed the fourth-order ERI tensor in terms of a set of third-order and second-order tensors of quadratically growing cardinality, where the number of elements of each tensor is small and independent of the system size. Compared to global DF, the complexity of computing the exact exchange regarding both CPU time and memory requirements is greatly reduced: First, the number of equations that have to be solved in the fitting procedure is equal to the number of negligible AO products,[183] which always scales linearly with system size. Second, the memory requirements are brought down to quadratic and further to linear when insignificant pair densities are not fitted. Finally, the exact exchange matrix is evaluated in N3: With eq , the elements of the exchange matrix K might be expressed aswhere P denotes elements of the density matrix. Cubic scaling is reached as it is possible to arrange all contractions in a way that never more than three atomic centers are involved. Taking into account distance effects, the scaling (with regard to both timing and memory) can be further brought down to linear as we will elaborate in section .

Numerical Considerations and Fit Set Quality

Compression of the ERI tensor via auxiliary basis set expansions usually introduces fitting errors ϵ in the ERIs (see (2)), where ϵ → 0 when the auxiliary basis approaches completeness. The rate of convergence of the auxiliary basis set expansion obviously depends strongly on the chosen fitting metric as well as on the nature of the fitting procedure: Robust fitting,[138]shows an error ϵ falling off bilinearly with the fitting error. However, due to its much lower computational complexity, we rely on nonrobust fitting instead,resulting in errors linear in the fitting error.[189] Consequently, due to the small number of auxiliary fit functions used to expand the pair densities (when compared to global fitting approaches), one would expect rather large errors in the computed integrals. However, the PARI-ERIs are usually very accurate approximations to the ERIs obtained without DF.[187] The exchange energy contribution (and consequently the HF/Hybrid-KS-DFT energy) from nonrobust fitting is unbounded from below,[183] making the SCF variationally unstable especially for large basis sets; calculations using quadruple-ζ-quality basis sets are often unreliable. The procedure becomes numerically more stable when rather large auxiliary fit sets are used.[186,189] In ADF, they are obtained as even-tempered series,[192−194] where the quality is controlled by the number of fit functions placed within a given range from the atomic center. Table shows the number of fit functions for representative atom types for the fit sets employed in this work.

Table 1

Number of Auxiliary Fit Functions (Angular Part Expressed in Real Spherical Harmonics) for Representative Types of Atoms for Different Fit Set Qualities on the All-Electron Level

	no. of auxiliary functions (composition)
quality	H	C	Au
Normal	62 (6s5p4d3f)	117 (10s9p5d4f3g)	779 (28s28p23d19f16g12h11i)
VeryGood	132 (10s7p6d5f4g)	209 (21s12p11d7f6g)	858 (31s22p20d19f16g16h16i)

DF-MP2 and AO-MP2 Equations

The canonical MP2 correlation energy EcorrMP2 for a closed shell molecule can be expressed aswhere (ia|jb) denotes an ERI in Mulliken notation, i, j (a, b) denote occupied (virtual) MO indices and ϵ (ϵ) denote diagonal elements of the occupied (virtual) blocks of the Fock matrix in the MO basis. Due to the transformation of the ERIs from the AO to the MO basis, the computational effort for the evaluation of this expression scales as N5. Using global DF, the prefactor of the computation of EcorrMP2 can be lowered considerably. The four-center integrals in the AO basis are then expressed aswith the third-order tensor C and second-order tensor V. With any local fitting approach, the MP2 energy can be evaluated in exactly the same way. Within the PARI approach one might simply obtain C and V from

In fact, employing PARI the prefactor of canonical DF-MP2 could be reduced further using sparse matrix algebra (the rank of the third-order tensor only scales linearly). In this work, we have not explored this possibility as far greater speed-up can be achieved working in the AO basis.

Applying the transformationsuggested by Häser and Almlöf[150−152] results inwithand

The sum over τ and the weight ωτ in (16)–(18) result from the evaluation of (15) by numerical quadrature. The optimal values for the set of quadrature points and weights can simply be obtained by least-squares minimization of the error distribution functionhowever, the minimax approximation[195,196] is a computationally more efficient approach. is obtained by transformation of the AO-ERI tensor according toemploying the PDMs P(τ) and Q(τ), given asThe quality of MP2 energies can generally be improved by empirically scaling individual contributions to it, giving rise to the popular spin-component-scaled (SCS)-[115,197−199] and spin-opposite-scaled-(SOS) MP2[116,117,200,201] approaches, where SCS-MP2 is often more accurate than SOS-MP2. However, for SOS-MP2 the second term on the right-hand side of (16) is completely neglected and the SOS-MP2 energy is obtained fromwhere usually cSOS = 1.3 is chosen.[116] This part of the MP2 energy can be evaluated with considerably lower computational cost than the same-spin part as it factorizes in a more favorable way and tensor contraction techniques can be used more efficiently.[164,180,181] Given the fact that DHs based on SOS-MP2 usually come very close to the accuracies of SCS-MP2 based ones,[127,202] the construction of fast SOS-MP2 methods alone seems to be highly desirable. Thus, in our efforts to develop a low-cost MP2 implementation, we have focused on the evaluation of only. The resulting algorithm and its implementation will be the subject of the next sections.

SOS-AO-PARI-MP2 Equations

In the AO basis, EcorrMP2 can be obtained from summing up the contributions from all pairs of atoms. Such a decomposition has already been suggested by Ayala and Scuseria nearly two decades ago;[203] here it arises quite naturally from the PARI approach. In particular, the Coulomb term can be obtained fromwhere denotes the contribution to e(2) from atom pair (A, B). Inserting (7) into (17) and dropping the index τ in the following, we obtain with the help of (20)Equation factorizes (unlike the corresponding expressions for ) according toand q (with q we will denote the set of all tensors qAA) is given asEquation can be computed in N3, and the same is true for (25): We first evaluate the local PDMs (21) for each pair of atoms. Then we half-transform the fit coefficients according to

Both transformations only involve three atomic centers, namely, A, B, and B′. Furthermore, f and g are two-center quantities and the memory required to store f and g (the set of all tensors f/g for all pairs of atoms (A, A′)) scales quadratically with the number of atoms. Subsequently, the first term in (25) is evaluated according toand the second term asresulting in

In fact, the chosen sequence of tensor contractions is closely related to a recent SOS-MP2 implementation by the Ochsenfeld group,[164] being even more apparent if one considers a monatomic system only for which the local pair fitting approach is equivalent to global DF. In all contraction steps, at most three atomic centers are involved, implying cubically scaling computation of the SOS-AO-MP2-PARI energy without any further consideration of distance effects. Furthermore, the three-center quantities f̃ and ẽ are evaluated on the fly, so that the memory requirements of the algorithm scale quadratically. As will be discussed in the next section, taking into account distance effects scaling of computation time and memory can be reduced further. We note that PARI also allows for quartic scaling computation of the exchange-like term in EcorrMP2. At this point we want to emphasize the strong similarity of our algorithm to the THC-approach by Martinez and co-workers:[174,175,180,204] Both methods exploit the locality in the atomic orbital basis set directly by decomposition of the ERI tensor into factor matrices that grow initially as N2, and as N as soon as a given system size is reached. Consequently, for the computation of EcorrMP2, the same formal scaling is reached: N3 for e(2) and N4 for e(2).

Distance Effects

As the basis functions χμ, χν are localized around their atomic centers, their overlap will decrease with increasing distance between the atoms on which they are centered. Consequently, the value of the overlap integral O (eq ) will approach zero with growing distance between the centers A and B. To exploit this behavior, we define a threshold (we will refer to this approximation as distant centers approximation for basis functions (DCAB) in the following) and consider a basis function as negligible for |r| > dμ ifwhere dμ is some basis function dependent effective radius to be determined at runtime. The procedure is illustrated in Figure .

Figure 1

Schematic illustration of the dependence of dμ on for two different types of functions. As a p-type function (blue) generally decays slower then an s-type one (red), its effective radius is larger.

Consequently, we only compute the fit coefficient c whenIn practice, we reorder the basis functions from the most diffuse (most slowly decaying) to the tightest one for each atom A, so that the dimension of the fit function tensor cAB approaches 0 forIfall tensor contractions involving the fit coefficient tensor corresponding to the pair (A, B) will be skipped. In the same way, we are also skipping tensor contractions involving the tensor VAB if the range of the Coulomb potential due to A does not overlap with any basis function on B (distant centers approximation for Coulomb potential (DCAC)). As the Coulomb potential decays only as |r|–1, this approach will only be effective for very large molecules. However, the Coulomb potential due to a pair density can be approximated using the well-known multipole expansion of the Coulomb operator.[205] Thus, the interaction between two pair densities ρ, ρ is evaluated via multipole expansion (recall that the fit functions are always assumed to be centered on the second atom of the pair) if

Similarly to the procedure for the basis functions, the actual values dα are controlled via a threshold . Considering multipole moments of up to lmax = 3, the dimension of the tensor VBD reduces to 16 × 16, for realistic fit sets being considerably lower than NfA × NfB (compare with Table ). The multipole approximation (MA) as well as the DCAB and DCAF approximation are used to speed up the computation of the post-SCF energy correction and the SCF itself. Exclusively for the MP2 part, we exploit sparsity of the density matrix in a way that we avoid contractions with half-transformed fit function tensors ifWe will refer to this approximation as neglect of half-transformed fit coefficient tensors (NHF) approximation. Clearly, as we do not use any localization techniques for the density matrix in our current implementation, this approximation can only be effective for spatially very extended systems. However, the density matrix and the half-transformed fit coefficient tensors show a high degree of sparsity, so that sparse matrix algebra techniques could efficiently be exploited here without a conceptual change of our implementation.

In practice, the efficiency of the possible screening options depends on thresholds, molecular geometry, and the diffuseness of the AO basis. Due to the rather large prefactor of Nf2Nb2, the scaling is dominated by steps 3a/3b, whereas the computational time for step 1 is always negligible. The asymptotic scaling of wall clock time and memory of our algorithm under consideration of screening effects is presented in Table .

Table 2

Outline of the Basic SOS-AO-PARI-MP2 Contraction Steps with Asymptotic Scaling (Big-O Notation Implied) and Memory Requirements under Consideration of Distanc Effectsa

step		asymptotic scaling	distance effects	memory
	calculate P and Q	N3		Nb,all2
2a	eν′μα = Pν′νcνμα
2b	fν′μα = Qν′νcνμα	N2	DCAB	Na2Nb2Nf (NaNb2Nf)
3a	mαα′ = eν′μαQμμ′cν′μ′α′
3b	nαα′ = fν′μαPμμ′cν′μ′α′	N	DCAB, NHF	on the fly
4a	gαα′ = eμ′μαfμμ′α′
4b	hα′α = fμμ′α′eμ′μα	N2		on the fly
5	qαα′ = mαα′ + nαα′ + gαα′ + (h)αα′T			Na2Nf2 (NaNf2)
6	Zαβ′ = qαα′Vα′β′	N2	MA, DCAF
7		N2

For each step the employed distance effects are given. The Einstein sum convention is used, which here always involves summation over the respective atomic centers. The memory requirements given in brackets refer to the memory saving variant of our algorithm.

A very important question regarding the feasibility of a local MP2 computation is its memory requirement. The half-transformed fit coefficient tensors f and g can be kept in memory for rather large molecules consisting of several hundreds of atoms. Although the memory requirements for these quantities is formally linearly scaling, they can hamper the application of our algorithm to very large systems. However, storage of f and g for each pair of atoms can be avoided if these quantities are recalculated prior to the contraction (28) (i.e., if step 2a/2b is repeated before 4a/4b) and with slight changes of the loop structure, storage of q can also be avoided. As the number of non-negligible fit coefficients and fit functions grows linearly, our algorithm is then order-N in memory. We also note that the practical memory bottleneck in our implementation is the storage of the untransformed fit coefficient tensors for compact systems. This can be attributed to our large auxiliary fit sets and we are planning further optimizations in this direction.

Loop Structure and Parallelization

Our algorithm is implemented by setting up two nested loops running over all pairs of atoms, which are closed whenever a quantity needs to be stored for each pair. Whenever we sum over a third center, a third loop over atomic shells is invoked. To be memory efficient, the loops need to be organized in a way that storage of three-index quantities is always avoided. The concept is demonstrated for the KIII contribution to the exact exchange matrix in Algorithm 1 (see eq ).

One might easily identify the first step in this algorithm as quadratically scaling, whereas the second one scales linearly. Although contractions involving V can rarely be skipped in practice, the multipole approximation considerably reduces the prefactor of this step. As the number of fit functions is usually larger than the number of basis functions, the second step determines the overall timing for the evaluation of KIII, explaining why we observe a subquadratic scaling behavior in practice (see section ). The SOS-MP2 correlation energy can be evaluated in a similar way. The algorithm for steps 2a/2b is outlined in Algorithm 2. The aforementioned quadratic asymptotic scaling of this step can be verified easily.

The loop structure of the presented algorithms suggests a parallelization strategy in which the tensor contractions associated with a certain pair of atoms (A, D) are distributed over all processes. Each tensor contraction is then performed on a single core. To avoid overhead for the MP2 part, we employ a two-level parallelization strategy where all numerical quadrature points are distributed over all available nodes and second the outermost two nested loops over atomic shells are distributed over all cores on this node. Thus, the number of nodes that can be used efficiently is limited by the number of numerical quadrature points Nnq (usually less than 10). However, we do not expect this to be a relevant issue for the most probable targets of our algorithm.

Benchmark Calculations and Discussion

Computational Details

We performed PARI-MP2 and SOS-AO-PARI-MP2 calculations with a locally modified development version of ADF.[137,206,207] We employed standard STO basis sets[190] and auxiliary fit sets for the evaluation of the exact exchange as described in section . If not indicated otherwise, all computations have been performed on the all-electron level.

To rule out the possibility of numerical issues throughout our benchmark calculations, the numerical integration quality, as well as the quality of the DF for the Coulomb part,[208] has been chosen to be better than default (Good quality, see refs (208) and (209)) if not stated otherwise. For the SCF, the mixed ADIIS+SDIIS method[210] has been employed.

All dimerization energies have been computed using the counterpoise (CP) method of Boys and Bernardi[211] to correct the basis set superposition error (BSSE). For all systems involving transition metals, relativistic effects have been treated with the zero-order regular approximation (ZORA)[212−215] in conjunction with ZORA-optimized basis sets and the minimum of neutral atomic potential approximation (MAPA).

If not stated otherwise, we used Nnq = 6 for all SOS-AO-PARI-MP2 calculations, where the numerical quadrature has been performed with a code[216] developed by Helmich-Paris.[195,196] For the evaluation of the exact exchange as well as for all SOS-AO-PARI-MP2 calculations, the default screening thresholds are

All calculations presented in this work were performed on a 2.2 GHz intel Xeon (E5-2650 v4) with 24 cores and 128 GB RAM. All binaries have been created using the GNU Fortran compiler.

Accuracy and Convergence with Basis Set Size

In this section we want to assess (a) the error introduced by the PARI-approach compared to that introduced by canonical DF-MP2, (b) the numbers of auxiliary fit functions required for accurate results as well as for a numerically stable SCF, and (c) the quality of our standard, non-correlation consistent STO basis sets. To this end, we performed benchmark calculations on different popular test sets: These are the s66 test set of weak intermolecular interactions[217] and test sets of relative conformational energies from the GMTKN30 database.[218] To assess the accuracy of our implementation in conjunction with an approximate treatment of relativistic effects, we calculated the HEAVY28 test sets of noncovalent interactions between heavy element hydrides.[219] We also report results for the L7[220] test set of weak intermolecular interactions with dimers of between roughly 50–120 atoms to asses the performance of our SOS-AO-PARI-MP2 implementation for large molecules. The entirety of all these test sets comprises of 187 data points.

We calculated the s66 test set of weak intermolecular interactions (London dispersion, hydrogen binding, π–π interactions) using TZP (triple-ζ with single shell of polarization functions) and TZ2P (triple-ζ with two shells of polarization functions) basis sets in conjunction with the Normal auxiliary fit set (see Table for the number of auxiliary fit functions for selected atoms). Figure shows deviations from our results to DF-MP2/CBS reference values[217] for all individual data points as well as mean absolute deviations (MAD) for each basis set. It is well-known that extrapolation to the basis set limit is only possible when correlation consistent basis sets of systematically increasing size are employed.[29] As a consequence, comparison of PARI-MP2/TZP(TZ2P) calculations to DF-MP2/CBS results cannot reveal if deviations can be attributed to basis set incompleteness or to the local pair fitting approach. To obtain a clearer idea about the reasons for the observed deviations from the DF-MP2/CBS results, we also compare our computed energies to DF-MP2/aug-cc-pVDZ and DF-MP2/cc-pVTZ reference values.[217] Both basis sets are comparable in size with TZP and TZ2P, with cc-pVTZ being the largest basis set with three shells of polarization functions, and TZP the smallest one, and the only one with only one shell of polarization functions. Figure clearly demonstrates that calculations on PARI-MP2/TZ2P/Normal-level yield results comparable to those of DF-MP2/aug-cc-pVDZ and DF-MP2/cc-pVTZ, with a MAD slightly better than DF-MP2/cc-pVTZ and slightly worse than DF-MP2/aug-cc-pVDZ. Also the sign-corrected maximum errors are with 2.8 kcal/mol for TZP and 2.3 kcal/mol for TZ2P in line with the two Gaussian type orbital (GTO) basis sets for which Figure also shows maximum errors considerably larger than 2 kcal/mol. Although the computed energies are on average very close to reproducing the DF-MP2/CBS references within a chemical accuracy of 1 kcal/mol, the PARI-MP2/TZP results are generally inferior to their TZ2P counterparts, which can safely be attributed to the smaller number of polarization functions.

Figure 2

Deviations from basis set extrapolated DF-MP2/CBS reference values[217] (in kcal/mol) for different basis sets for each data point in the s66[217] test set of weak intermolecular interactions. The aug-cc-pVDZ and cc-pVTZ reference values[217] have been computed with DF-MP2, whereas the TZP and TZ2P values have been obtained with our PARI-MP2 code, using the Normal fit set for both HF and MP2. MADs are given with respect to DF-MP2/CBS.

From the good agreement of our PARI-MP2 results with DF-MP2 for basis sets of comparable size, we conclude that the PARI-approximation does not seem to considerably degrade the accuracy of DF-MP2 dimerization energies. Furthermore, our findings suggest that our non-correlation consistent STO-type basis sets can compete with larger correlation consistent GTO-type basis sets with more shells of polarization functions.

A factor, not having been discussed so far, is the size of the auxiliary fit sets. As already pointed out, a lager fit set should help to ensure variational stability of the HF energy in the SCF. To this end, we computed the relative conformational energies in the GNTKM30 database (the ACONF, SCONF, PCONF, CyCONF, and ISO34 subsets) with Normal and VeryGood fit set for both HF and PARI-MP2 and also recalculated the s66 dimerization energies with the larger auxiliary fit set. Somewhat surprisingly, we found that the larger auxiliary fit set sometimes led to deteriorated PARI-MP2 energies, especially for the S66 test set, whereas the HF energies remained essentially unchanged (as expected). This observation is directly opposed to the fact that larger fit sets should generally improve any DF approximation.

We then carried out calculations where we retained the VeryGood fit set for the SCF but used the Normal fit set for the evaluation of the MP2 energy correction; i.e., we recalculated the ERIs after convergence of the SCF using a smaller fit set. The results are summarized in Figure .

Figure 3

PARI-MP2 results (in kcal/mol) for selected test sets of isomerization energies from the GNTKM30 benchmark set[218] as well as for the s66 test set. The MADs for each data set with respect to DF-MP2/CBS, as well as the MADs for the entirety of all test sets (in total 152 data points) are given. Key: basis set/HF auxiliary fit set//MP2 auxiliary fit set.

Clearly, for almost all test sets, the computed energies are nearly independent of the auxiliary fit set employed in the SCF. This is reflected in the MADs over all test sets, where literally no difference can be observed between Normal and VeryGood. Only for the SCONF test set, the VeryGood fit set is slightly inferior to Normal.

As recalculation of the ERIs after the SCF seems to cure the problems in the computation of the PARI-MP2 energy for the VeryGood fit set, one must conclude that neither the MOs nor the orbital energies resulting from the use of this fit set are problematic for the computation of EcorrMP2, but rather the fitting errors in the ERIs themselves. Due to the fact that for some test sets the VeryGood fit set yielded (rather small) improvements over the Normal one, we do not suspect a fundamental issue with the PARI approximation here but rather a numerical one. Clearly, the large number of auxiliary fit functions for the expansion of each pair of basis functions leads to linear dependencies in the fit set that might cause numerical problems due to overfitting. We emphasize that the VeryGood fit set is usually not inferior to the Normal one. However, the latter one seems to be numerically more stable for PARI-MP2 calculations and also completely sufficient for the purpose of the present study. Optimizing fit sets specifically for correlation methods seems to be highly promising for even more accurate PARI-MP2 energies but is out of the scope of this work.

The clearly improved results over the TZP basis set when TZ2P is used for the s66 test set can generally not be observed for the test sets of conformers. For the ACONF (alkane conformers) and PCONF (tripeptide conformers) test sets, the TZP energies are even slightly better than the TZ2P ones. This can only be explained with error cancellation: It is well-known that MP2 often tends to overestimate correlation energies,[221] so that a more incomplete basis set might yield better results for certain systems. Only for CYCONF (cysteine conformers) and ISO34 (isomerization energies of organic molecules), a clear improvement over TZP can be observed when TZ2P is used instead. As the overall MADs over all 152 data points reveal, all considered combinations of basis and fit sets reproduce the DF-MP2/CBS reference within chemical accuracy on average, where the TZ2P basis set is in many (but not all) cases superior to TZP.

As shown in Table , we find only moderate maximum errors for the CYCONF and ACONF test sets with 0.5 and 0.7 kcal/mol, respectively for the TZP basis set and 0.3 and 0.9 kcal/mol, respectively, for the TZ2P one. For the SCONF and PCONF test sets we find larger maximum errors of over three kcal/mol irrespective of basis set and fit set and nearly 4 kcal/mol for the ISO34 test set with the TZP basis set, reflecting the slow basis set convergence of MP2 correlation energies.

Table 3

Maximum Sign Corrected Errors for the S66 Test Set and the Test Sets of Relative Conformational Energies for the Herein Investigated Combinations of Basis Sets and Fit Sets (All Energies in kcal/mol)a

test set/fit set	ACONF	CYCONF	ISO34	SCONF	PCONF	S66
TZP/VeryGood//Normal	0.68	0.50	3.86	3.20	3.46	2.83
TZ2P/VeryGood//Normal	0.93	0.28	3.62	3.09	3.79	2.13
TZP/Normal	0.67	0.53	3.87	3.20	3.45	2.87
TZ2P/Normal	0.92	0.25	3.56	3.08	3.76	2.25

First column: basis set/HF auxiliary fit set (//MP2 auxiliary fit set).

As our findings show that the more expensive VeryGood fit set yields no improvement over the Normal one, we strongly recommend using the latter one. Considering error cancellation between e(2) and e(2) as highly unlikely, our findings also apply to all spin-scaled variants of MP2.

For heavy elements we found excellent agreements to the CCSD(T)/CBS reference values for the HEAVY28 test set.[219] Already on the TZP/Normal/ZORA-level, our results show a MAD of only 0.30 kcal/mol. This is actually better than DF-MP2/CBS for which a MAD of 0.41 kcal/mol has been reported when small effective core potentials of the Stuttgart/Cologne group[222,223] for elements with Z > 36 are used.

After having assessed the performance of our MO based PARI-MP2 algorithm, we finally explore the accuracy of our SOS-AO-PARI-MP2 implementation for large molecules. To this end, we computed the dimerization energies from the L7 test set with varying thresholds controlling the distance effects outlined in section . As we are not aware of SOS-MP2 reference energies for this test set, we compare our results to the available QCISD(T)/CBS values.[217] We emphasize that the only approximation in SOS-AO-PARI-MP2 that is not already present in SOS-PARI-MP2 is the numerical quadrature (15). For an account on its convergence with respect to Nnq we refer to the literature.[164] We only note that Nnq = 6 is usually sufficient to achieve millihartree accuracy. However, a broader range of orbital energies and a smaller HOMO–LUMO gap require more quadrature points and Nnq = 8 might be more appropriate. To rule out the possibility of inaccuracies due to a too small Nnq, we have decided to use Nnq = 8 for the numerical computation of (15).

First, we tested the performance of our method for different thresholds, controlling the distance effects described in section . We do not consider variations of ϑNHF here and also found it useful to group the remaining thresholds in tiers. For each of these tiers, as described in Table , we have calculated the L7 test set with the TZP basis set and the Normal auxiliary fit set. Our results are shown in Figure , where the upper part shows the deviations of computed reaction energies with respect to the Basic tier. Going up to the Normal tier (painted area), the largest change in energies occurs for c3gc with roughly 0.3 kcal/mol. Proceeding from Normal to Good, the energy does only change marginally for all dimers, justifying the Normal tier as our default. Essentially no difference can be observed between the last two tiers.

Table 4

Definition of Threshold Tiers

	Basic	Normal	Good	VeryGood
ϑDCAB	2 × 10–3	1 × 10–3	3 × 10–4	1 × 10–4
ϑDCAC	2 × 10–2	1 × 10–2	1 × 10–3	1 × 10–4
ϑMA	3 × 10–2	3 × 10–3	3 × 10–3	3 × 10–4

Figure 4

Upper part: Deviations from results obtained with the Basic tier of thresholds for the individual reaction energies. Lower panel: MADs with respect to the QCISD/CBS reference[217] for SOS-AO-PARI-MP2 calculations, as well as for DF-MP2/CBS. All energies are in kcal/mol. The naming of the dimers follows Řezáč et al.[217]

As shown in the lower half of Figure , using smaller thresholds leads to an artificial improvement of the dimerization energies: MP2 often exaggerates correlation effects as its double excitations do not couple.[221] Lowering the cutoff thresholds, however, corresponds to the neglect of correlation from distant pairs. Thus, the underestimation of dispersion interactions due to the neglect of long-range correlation effects partly compensates for the overestimation of the correlation energy within MP2, explaining the trend in the observed reaction energies.

We finally note that the TZ2P basis set gives improvements over TZP, although the MAD is still higher than 4 kcal/mol. However, this is significantly better than full DF-MP2/CBS calculations for which a MAD of 6.58 kcal/mol has been reported in the literature.[217]

To conclude this section, we think that the presented data clearly demonstrates (i) that our default fit set is completely sufficient to compute accurate PARI-MP2 energies in a numerically stable way, (ii) that the deviation from the DF-MP2/CBS reference values can (at least to a great extent) be attributed to the basis set error, (iii) that our PARI-MP2 implementation used with non-correlation consistent basis sets of triple-ζ quality yields errors comparable to the ones from DF-MP2 with correlation consistent basis sets of the same size, and finally, (iv) that our implementation also yields accurate and reliable energies for large systems when rather conservative default screening thresholds are used.

Performance and Timing

We analyzed the performance of our SOS-AO-PARI-MP2 implementation on a series of linear alkane chains as an optimum-case for local correlation methods. Our results are summarized in Table .

Table 5

CPU Times and Scaling Behavior with Respect to the Systems Size Relative to the Previous Calculation (in Parentheses) in Terms of the Polynomial Coefficient x in N for SOS-AO-PARI-MP2 Calculations on Linear Alkane Systems Using TZP and TZ2P, Respectively, and Normal Fit Set Quality (All Calculations on a Single Core, All Timings in min)

		timing
n (CnH2n+2)	no. of bf	total		MP2 alone		MP2 time [% of full calc]
TZP
20	632	23.5		3.9		16.6
40	1252	76.8	(1.70)	17.5	(2.16)	22.8
80	2492	255.3	(1.79)	75.3	(2.11)	29.5
160	4972	907.3	(1.83)	323.6	(2.10)	35.6
TZ2P
20	1084	41.7		8.4		20.1
40	2144	139.4	(1.74)	37.4	(2.15)	26.8
80	4264	473.0	(1.76)	156.0	(2.06)	33.0
160a	8504	1811.2	(1.94)	752.3	(2.27)	41.5

Calculated with the more memory efficient variant of the algorithm.

In all computations, the evaluation of the MP2 correlation energy only accounts for between 20 and 42% of the total wall clock time. Thus, the overall scaling is dominated by the HF part, reaching subquadratic scaling already for the shortest chains considered here, whereas for the computation of the post-SCF energy correction quadratic scaling is observed for both basis sets. This can be attributed to the rather conservatively chosen screening thresholds due to which NHF screening becomes practically irrelevant, even for large and spatially extended systems. It should also be emphasized that the scaling is not strongly affected by the number of diffuse functions in the basis set. For the TZ2P basis set containing a larger number of diffuse polarization functions, the scaling is even slightly better than for the TZP basis set. Due to the huge memory requirements to store f and g, we have been forced to switch to the slightly slower, but more memory efficient variant of our algorithm for the computation of C1600H322, explaining the comparatively high increase in wall clock time for this step.

The excellent scaling behavior of our method for relatively small systems is also reflected in the early crossover point with respect to canonical SOS-PARI-MP2; for C10H22, the SOS-AO-PARI-MP2 energy correction alone is computed in 45 s on the TZP-level of theory, while the respective canonical calculation already takes 93 s.

On the example of C40H82 in the TZP basis, we give an estimate on the efficiency of our parallelization strategy: Comparing wall clock times obtained with 1, 8, and 24 cores on the same node, we find parallel speedups of 6.4 and 15.8, respectively. Parallel timing results are also presented for stacks of backbone-free DNA in Table .

Table 6

CPU Times and Scaling Behavior with Respect to the Systems Size Relative to the Previous Calculation (in Parentheses) in Terms of the Polynomial Coefficient x in N for SOS-AO-PARI-MP2 Calculations on Backbone-Free DNA Stacks on the TZP/Normal Level of Theory (All Calculations on a Single Node with 24 Cores, All Timings in min)

no. of units	no. of bf	total		MP2 alone		MP2 time [% of full calc]
1	848	4.1		0.7		17.1
2	1696	20.1	(2.29)	4.3	(2.62)	21.2
4	3392	79.7	(1.99)	20.9	(2.28)	26.2
6	5088	185.1	(2.08)	52.1	(2.25)	28.1

Here, each unit consists of an adenine–guanine and a cytosine–thymine base pair, separated by 3.4 Å. Although these systems are still spatially extended, they are considerably more compact as linear alkanes. For the largest of these systems considered here with 354 atoms and 5088 basis functions, the SOS-AO-PARI-MP2 energy can be evaluated in approximately 3 h on a single compute node, with the computation of the SOS-MP2 energy only accounting for 28% of the total elapsed time. Although the overall scaling is slightly worse than for the linear alkane chains, one still discovers the onset of subquadratic scaling for the whole calculation and convergence to quadratic scaling for the SOS-MP2 calculation only.

The efficiency of our implementation is illustrated on different types of realistic, compact systems (see Figure ) in Table , where we give detailed timings for the most expensive steps of the SOS-AO-PARI-MP2 calculations from Table .

Figure 5

Realistic 3D systems employed in this work. Upper panel from left to right: 4b, 7b from the S30L testset of Grimme and co-workers,[224] and a (H2O)142 water cluster[156] from the Ochsenfeld benchmark set. Lower panel from left to right: A (S8)20 sulfur cluster,[156] a DNA segment from adenine–thymine base pairs[159] (both from the Ochsenfeld benchmark set), and a substituted cluster of 21 Au atoms from Jones et al.[225]

Table 7

Detailed Wall Clock Times (in min) for SOS-AO-PARI-MP2 Calculations on Selected Realistic 3D Systems on the TZP/Normal Level of Theory on a Single Node with 24 Coresa

	4b	7b	(H2O)142	DNA4	(S8)20	Au21S(SCH3)15b,c
no. of atoms	158	153	426	260	160	97
no. of bf	2768	2248	4544	3638	4480	2414
Timings
total	76.2	40.3	186.0	104.3	177.9	269.2
total MP2	15.6	9.4	56.7	28.7	75.5	106.7
step 1d	0.06	0.04	0.21	0.10	0.19	0.05
step 2a/2b	2.4	1.4	6.63	3.6	11.8	19.62
step 3a/3b	10.0	5.7	36.0	16.4	49.6	75.6
step 6	2.2	1.6	18.0	6.8	6.2	6.3

The first two structures are taken from the S30L test set,[224] structures 3–5 are from the test set of Ochsenfeld and co-workers[156,159] and the structure of the last molecule has been taken from Jones et al.[225]

Relativistic effects have been treated on the ZORA/MAPA level of theory.

Due to the small HOMO–LUMO gap, Nnq = 8 was chosen.

Numbering of steps refers to Table .

The most expensive part for each calculation is the SCF. The wall clock time for the calculation of the SOS-AO-PARI-MP2 energy correction is clearly dominated by step 3a/3b, whereas the cubic scaling evaluation of the PDMs (step 1) is negligible. Step 3a/3b is also the part of our algorithm that would probably profit most from the exploitation of sparsity in the half-transformed fit coefficient tensors. For the very compact systems (S8)20 (compund e in Figure ) and Au21S(SCH3)15 (f) in Figure ), the half-transformation of the fit coefficients (step 2a/2b) also consumes a considerable share of the total wall clock time as distance effects do not come into play here. As the majority of the individual tensor contractions scale as Nb2Nf2 in the SCF and in the MP2 part, the large number of basis functions and auxiliary fit functions for the gold atoms (see Table ) makes the computation for this system particularly slow. The efficiency of step 6 is more or less independent from the molecular geometry, best seen on the example of the water cluster, indicating that the multipole expansion does not lead to large computational savings for these compact systems.

In the end, we comment on the memory requirements of our code. The computation of the SOS-AO-PARI-MP2 energy for C160H322 in TZ2P quality with roughly 8500 AOs could only been achieved by using our more memory efficient implementation, even when 128 GB memory is used and calculations on even larger systems are impossible due to the memory requirements of the SCF. Although we do not think that molecules of this size will be the main target of our implementation, this is currently a severe drawback of our algorithm and we are planning further optimization in this direction.

Conclusion

We have demonstrated on test sets of collectively 187 data points that dimerization energies and conformational energies from PARI-MP2 deviate from their DF-MP2/CBS counterparts by less than 1 kcal/mol on average, when non-correlation consistent STO-type basis sets of triple-ζ quality and our default auxiliary fit sets are used. We have also demonstrated the accuracy of this approach for large systems of more than 100 atoms and shown that our implementation reproduces CCSD(T)/CBS reference values better than DF-MP2/CBS for the HEAVY28 test set of noncovalent interaction energies between heavy element hydrides when relativistic effects are approximated on the ZORA/MAPA level.

Comparison to DF-MP2 calculations on the S66 test set shows that the error of our calculations is of the same order of magnitude as the basis set incompleteness error of GTO-type basis sets of comparable size. The maximum deviation observed is below 4 kcal/mol for the TZP basis set with only a single polarization function per atom and considerably lower than 1 kcal/mol for some of the investigated test sets. We expect significant improvement of these values by employing fit sets optimized for correlation methods as it is common practice for DF-MP2 with GTOs.[142,226−228] To the best of our knowledge, such fit sets have not been designed for STO based PARI yet and research in this direction is currently being pursued by our group.

Additionally, we have presented a quadratic scaling SOS-AO-PARI-MP2 algorithm. The overall evaluation of the SOS-MP2 energy scales quadratically and the post-SCF energy correction is computed considerably faster than the SCF itself for all system considered herein. Among others, we have demonstrated the efficiency of our approach on a very compact cluster of 160 sulfur atoms with 4480 basis functions, and a cluster of 142 water molecules with 4544 basis functions: Each all-electron calculation could be performed in approximately 3 h on a single compute node. Another attractive feature of our implementation is its early crossover to canonical SOS-PARI-MP2 in terms of wall clock time: For a linear alkane chain with only 10 carbon atoms, our AO based algorithm is already twice as fast as our MO based one. As a consequence, our algorithm is fast for medium and large compact systems of up to several hundreds of atoms, the bottleneck, with regard to both memory and wall clock time, being the SCF rather than the MP2 calculation. As we are mostly avoiding disk I/O, our approach is particularly appealing if the calculations are run on a machine with relatively slow disks, such as the nowadays ubiquitous computer clusters with external storage devices.

At the moment, the post-SCF energy correction does not scale as favorably as the SCF itself. Using sparse matrix algebra, we could possibly turn our current implementation into a truly linear scaling one. However, one might legitimately argue that this would not bring much additional value. Even if one could speed up the MP2 calculation alone by a factor of 2, the overall wall clock time for the computation of the SOS-MP2 energy for the sulfur cluster (where the MP2 is particularly slow compared to the SCF) would only be reduced by 20%. In other words, major efficiency gains can only be achieved when the bottleneck of the computation, the computation of the exact exchange in the SCF, is optimized.

Although larger computational savings could possibly be achieved for systems much larger than the ones presented herein, we do not think that even then these large systems would be a probable target of our algorithm. Furthermore, application of our algorithm to these systems is hampered by its rather unfavorable memory requirements (if one wants to avoid disk I/O) and further development of the code will rather focus on improvements in this direction. We are also planning to extend our implementation to periodic systems and to implement gradients as well.

Although we think that accurate and fast MP2 correlation energies are highly desirable in themselves, they are arguably most relevant in the framework of double-hybrid density functional approximations. At the moment, we are working on a comprehensive benchmark of state-of-the-art SOS based double hybrids for small as well as for large molecules. As shown only recently,[127] this class of double hybrids is in almost all cases not inferior to SCS-MP2 based ones, especially if the recently developed D4 extension[229,230] to the D3 dispersion correction[219,231,232] is used.