Theoretical Chemistry, Vrije Universiteit, De Boelelaan 1083, NL-1081 HV,
Amsterdam, The Netherlands Theoretical Chemistry, Vrije Universiteit, De
Boelelaan 1083, NL-1081 HV, Amsterdam, The Netherlands

# Towards Pair Atomic Density Fitting for Correlation Energies with Benchmark
Accuracy

Edoardo Spadetto<EMAIL_ADDRESS>Software for Chemistry and Materials NV,
NL, 1081HV, Amsterdam, The Netherlands  Pier Herman Theodoor Philipsen
<EMAIL_ADDRESS>Software for Chemistry and Materials NV, NL, 1081HV,
Amsterdam, The Netherlands  Arno Förster<EMAIL_ADDRESS>Software for
Chemistry and Materials NV, NL, 1081HV, Amsterdam, The Netherlands  Lucas
Visscher

###### Abstract

Pair atomic density fitting (PADF) has been identified as a promising strategy
to reduce the scaling with system size of quantum chemical methods for the
calculation of the correlation energy like the direct random phase
approximation (RPA) or second-order Møller-Plesset perturbation theory (MP2).
PADF can however introduce large errors in correlation energies as the two-
electron interaction energy is not guaranteed to be bounded from below. This
issue can be partially alleviated by using very large fit sets, but this comes
at the price of reduced efficiency and having to deal with near-linear
dependencies in the fit set. One posibility is to use global density fitting
(DF), but in this work, we introduce an alternative methodology to overcome
this problem that preserves the intrinsically favourable scaling of PADF. We
first regularize the Fock matrix by projecting out parts of the basis set
which gives rise to orbital products that are hard to describe by PADF. After
having thus obtained a reliable self-consistent field solution, we then also
apply this projector to the orbital coefficient matrix to improve the
precision of PADF-MP2 and PADF-RPA. We systematically assess the accuracy of
this new approach in a numerical atomic orbital framework using Slater Type
Orbitals (STO) and correlation consistent Gaussian type basis sets up to
quintuple-$\zeta$ quality for systems with more than 200 atoms. For the small
and medium systems in the S66 database we show the maximum deviation of PADF-
MP2 and PADF-RPA relative correlation energies to DF-MP2 and DF-RPA reference
results to be 0.07 and 0.14 kcal/mol respectively. When the new projector
method is used, the errors only slightly increase for large molecules and also
when moderately sized fit sets are used the resulting errors are well under
control. Finally, we demonstrate the computational efficiency of our algorithm
by calculating the interaction energies of large, non-covalently bound
complexes with more than 1000 atoms and 20000 atomic orbitals at the
RPA@PBE/CC-pVTZ level of theory.

###### keywords:

RPA, MP2, algorithms, numerical atomic orbitals, non-covalent interactions

##  1 Introduction

There is great scientific and commercial interest to successfully predict the
electronic structure of molecules and materials. Towards this aim, Density
Functional Theory (DFT)1, 2, 3 and the Hartree–Fock (HF) method 4, 5, 6, 7, 8
are indispensable tools, but often they capture exchange-correlation effects
only insufficiently. For instance, dispersion and polarization effects which
derive mainly from mid- to long-range electron correlation9 are not accounted
for by standard DFT functionals.10, 11. While empirical dispersion corrections
have been highly successful in describing some of these aspects of electron
correlation,12, 9 methods based on coupled cluster (CC)13, 14, 15, 16, 17
theory are generally considered to be the most precise class of correlation
methods and have been the workhorse for many high precision calculations.18
However, given their large computational cost,19 their usage is in practice
still limited to relatively small systems. Using massively parallel
implementations20, 21 and/or local approximations22, 23, 24, 25 it is in
principle possible to treat larger molecules.26, 27 Massively parallel
calculations do, however, require large computational resources, while the
errors introduced by the local approximations can give rise to uncertainties
of the order of several $\nicefrac{{kcal}}{{mol}}$ in some cases.28, 29, 26

An attractive alternative is provided by double hybrid (DH) density
functionals30, 31 which usually offer a good compromise between accuracy and
computational effort.32, 33, 34, 35 DHs combine DFT with methods which treat
correlation effects explicitly, mostly second-order Møller-Plesset
perturbation theory (MP2)36 but (more recently)37, 38, 39 also the random-
phase approximation (RPA) 40, 41 has been considered for this task. While MP2
is not necessarily very accurate and limited in its applicability, RPA is
gaining popularity42, 43, 44, 45, 46, 47, 48, 49, 50. RPA is applicable to a
wider class of systems than MP2, as it is unaffected by some of the drawbacks
of MP2, namely divergences for metals, small band gap systems40, 51 and large
molecules.52. Without increasing the computational effort one can also greatly
improve upon the accuracy of RPA by using $\sigma$-functionals.53, 54, 55
Another popular alternative is the inclusion of screened exchange which comes
essentially at MP2 cost.56, 57, 58, 59

RPA and MP2 are typically implemented with $\mathcal{O}\left(N^{4}\right)$60
and $\mathcal{O}\left(N^{5}\right)$ scaling with system size using global
density fitting (DF)61, 62, 63, 64 (DF-RPA60 and DF-MP2,64, 65 respectively).
Efficient implementations of these methods enable applications to systems with
a few hundred of atoms even at the quadruple-$\zeta$ level,50 but larger
systems are out of reach on standard hardware. For this reason, more efficient
algorithms and approximate implementations have been developed to improve the
scaling of both RPA and MP2. Common strategies are the usage of localized
orbitals,66, 67 cluster-in-molecule (CIM) approaches,68, 69 or implementations
which rely on sparsity in the atomic orbital basis.70, 71, 72, 73, 74, 75, 76,
77, 78, 79, 80, 81, 82, 83, 84, 85, 86 In the latter class of methods,
implementations using local DF approximations have gained increasing
popularity.76, 82, 83 While they do not achieve linear scaling with systems
size, they typically come with a very small prefactor83 and are believed to
only introduce minor errors compared to canonical, molecular orbital based
implementations.83, 87

One particular flavour of local DF approximations is Pair Atomic Density
Fitting (PADF),62, 88, 89, 90, 83 also known also as pair atomic resolution of
the identity (PARI),91, 83 concentric atomic density fitting,92, 93 or RI-
LVL.94 PADF has originally been introduced to speed up the construction of the
Hartree contribution in non-hybrid DFT calculations but was later generalized
to accomodate also the formation of the exact exchange matrix in HF or hybrid
DFT calculations. For a comparison of PADF-HF to other approximate exact
exchange algorithms see ref. 95. PADF can also be used to reduce the
asymptotic scaling of RPA and Spin-opposite scaled (SOS)-MP296 to formally
cubic.83 However, quadratic scaling is often observed in practice since the
prefactor of the cubic steps is small.83

This speedup comes at the cost of errors which can in principle cause
variational collapse of HF calculations to solutions corresponding to
artificially low energies.90 In practice, this is usually only an issue when
insufficiently large fit sets are employed. What often is more problematic
that PADF also leads to an artificial increase of the magnitude of correlation
energies. Unless unrealistically large fit sets are used, this is difficult to
avoid these errors and this can then also affect the precision of relative
energies that are the typical target of quantum chemical calculations.

In the Amsterdam modelling suite (AMS)97, 98, the issues of PADF are mitigated
by applying a projector to the exact HF exchange matrix in order to prevent
variational collapse. The same projector can also be applied to orbital
coefficients in order to reduce errors in post SCF methods. This strategy has
already been successfully employed in the past for many-body perturbation
theory (MBPT) based calculations.59, 99 However, systematic benchmarks against
other codes using the same basis sets were not yet performed. For this
purpose, we report here an implementation of PADF-MP2 and PADF-RPA in the
numerical atomic orbital (NAO) based code BAND.100, 101, 102, 103 This
implementation allows us to use Gaussian type orbitals (GTO) as basis sets and
therefore to systematically investigate the accuracy of the PADF-MP2 and PADF-
RPA implementations in AMS for relative correlation energies with respect to
global DF based implementations (DF-MP2, DF-RPA) in Psi4 and TURBOMOLE 104,
105. Similar benchmarks of PADF based correlation energies have already been
reported by Ihrig et al. using the FHI-AIMS code106, 107, 108, 109 who found
excellent agreement of PADF-MP2 and PADF-RPA with DF-MP2 and DF-RPA94 and also
by Tew.110 However, these authors focused on small and medium sized molecules
only. To assess whether this accuracy also holds for larger systems and large
basis sets, we herein report benchmarks for non-covalently bound dimers with
up to 200 atoms and for large GTO-type basis sets up to quintuple-$\zeta$
(5Z). In our benchmarks, we focus exclusively on non-covalent interactions.
This is mostly due to the availability of accurate reference values for
datasets containing large molecules, like the L7111 or the S30L112
compilations.

This paper is organized as follows: In section 2 we review the PADF method and
introduce the projector method (PM). We also sketch how PADF can be used to
achieve low-order scaling implementations of RPA and SOS-MP2. For more
details, we refer to previous work.83, 113 After an outline of our
computational details in section 3, we assess the accuracy of relative PADF-
MP2 and PADF-RPA correlation energies in section 4. Our calculations show that
PADF-SOS-MP2 is in excellent agreement to DF-SOS-MP2. Our interaction energies
for the S66 dataset114 show maximum absolute deviations for PADF-MP2 and PADF-
RPA with respect to the reference results of 0.07 and 0.14 kcal/mol
respectively, irrespective of the chosen basis set. For much larger molecules,
we observe only a negligible loss in accuracy for SOS-MP2 and we find the PM
to be decisive to obtain accurate results. The loss in accuracy is more
pronounced for PADF-RPA, but errors are smaller than errors due to basis set
incompleteness or due to local correlation approximations for large
molecules27, 26 To showcase the efficiency of our implementation, we calculate
PADF-RPA interaction energies of eight large non-covalently bound complexes at
the triple-$\zeta$ (TZ) level, with up to 1000 atoms and more than 20000 AOs.
Finally, section 5 summarizes and concludes this work.

##  2 Theory

Throughout this paper, the indicies $\\{i,j,...\\}$ ($\\{a,b,...\\}$) refer to
occupied (virtual) orbitals, and the indices $\\{p,q,...\\}$ refer to general
molecular orbitals. Primary basis functions are labeled with
$\\{\mu,\nu,\kappa,\lambda,...\\}$ while $\\{\alpha,\beta,\gamma,...\\}$
denote fit functions. $\\{A,B,C,...\\}$ denote atoms. $o$ is a generic index
which can either denote a primary basis function or a fit function.

###  2.1 Density Fitting

We use Mulliken notation throughout this work, in which the generic form of
two-electron integrals is given by

$\mathcal{K}_{\mu\nu\kappa\lambda}=\int
d\bm{r}d\bm{r}^{\prime}\chi^{*}_{\mu}(\bm{r})\chi_{\nu}(\bm{r})\mathcal{K}(\bm{r},\bm{r}^{\prime})\chi^{*}_{\kappa}(\bm{r}^{\prime})\chi_{\lambda}(\bm{r}^{\prime})\;,$
(1)

with $\mathcal{K}(\bm{r},\bm{r}^{\prime})$ being a general (non-local) kernel.
Important examples for $\mathcal{K}$ are the electron-electron interaction,
$v_{c}(\bm{r},\bm{r}^{\prime})=\frac{1}{|\bm{r}-\bm{r}^{\prime}|}$ which is a
key ingredient in HF and post-HF methods, and the non-interacting
polarizability $P(\bm{r},\bm{r}^{\prime})$ (for a certain value of imaginary
frequency or time) which appears for instance in RPA or SOS-MP2. The symbol
$\chi_{\mu}$ refers to an atom-centered basis function which belongs to a
basis set of $N_{\text{bas}}$ functions $\\{\chi_{\mu}(\textbf{r})\in X\hskip
4.0pt\forall\hskip 4.0pt1\leq\mu\leq N_{\text{bas}}\\}$. We assume these
functions to be composed of an angular part $Y_{l}^{m}(\theta,\phi)$ and a
radial function $R_{n}(|\textbf{r}_{A}|)$,

$\chi_{\mu}(\bm{r})=\chi_{lmn,A}(\theta,\phi,|\textbf{r}_{A}|)=Y_{l}^{m}(\theta,\phi)R_{n}(|\textbf{r}_{A}|)\;.$
(2)

The radial part only depends on the distance from atom $A$, $\textbf{r}_{A}$.
The angular part $Y_{l}^{m}$ is a spherical harmonic function with angles
defined in the local coordinate system of atom A.

Representing $\mathcal{K}$ in this way, the memory required to store all the
integrals defined by (1) grows as $\mathcal{O}\left(N^{4}\right)$ with system
size. Furthermore, evaluating these electron-electron interaction integral
explicitly is difficult, when sets of STOs or NAOs are chosen as primary
basis. It is therefore convenient to look in more detail at the set of
functions $F_{p}=\\{\chi_{\mu}(\bm{r})\chi_{\nu}(\bm{r})\hskip
4.0pt\forall\hskip 4.0pt1\leq\mu\leq N_{\text{bas}}\text{ and }1\leq\nu\leq
N_{\text{bas}}\\}$ that is obtained by gathering all unique products of two
basis functions and investigate whether this function set can be represented
in a more compact form via density fitting.

We first define $\hat{\mathcal{K}}$ formally115 as a linear operator in a
Hilbert space $\mathcal{H^{K}}$: $\hat{\mathcal{K}}f(\bm{r})=\int
d\bm{r}^{\prime}\mathcal{K}(\bm{r},\bm{r}^{\prime})f(\bm{r}^{\prime})$, with
the inner product on $\mathcal{H^{K}}$ defined as

$\displaystyle(f|g)=\int f(\bm{r})\;\hat{\mathcal{K}}\;g(\bm{r})d\bm{r}.$ (3)

The set of fit functions that is defined in (PA)DF forms the basis $F_{f}$ and
spans a subspace $X_{f}$ of the full Hilbert space $\mathcal{H^{K}}$. Given
$f_{\alpha}(\bm{r})$ and $f_{\beta}(\bm{r})\in F_{f}$, the expansion of
$\hat{\mathcal{K}}$ in $F_{f}$ is

$\mathcal{K}_{\alpha\beta}=(f_{\alpha}|f_{\beta}).$ (4)

Likewise we have that the basis $F_{p}$ spans the space
$X_{p}\subset\mathcal{H^{K}}$, with the integrals taking the form (1) or more
concisely $(\chi_{\kappa}\chi_{\lambda}|\chi_{\mu}\chi_{\nu})$. If
$X_{p}\subseteq X_{f}$, $\chi_{\mu}(\bm{r})\chi_{\nu}(\bm{r})$ can be
expressed exactly in terms of fit functions and we can consequently use the
compact expression (4) instead of (1). In practice we find that part of the
product space is not spanned by the fit functions ($X_{p}\setminus
X_{f}\neq\emptyset$). To characterize errors made by fitting basis functions
products with the fit set $F_{f}$ we therefore write members of the product
basis $F_{p}$ as

$\displaystyle\chi_{\mu}(\bm{r})\chi_{\nu}(\bm{r})=\sum_{\alpha}f_{\alpha}(\bm{r})c_{\mu\nu\alpha}+e_{\mu\nu}(\bm{r})$
(5)

where $e_{\mu\nu}(\bm{r})$ is an error function which accounts for the fact
that $X_{f}$ does not completely span $X_{p}$. To keep the notation short we
only indicate explicitly the dependence of this error function of the basis
function pair indices $\mu$ and $\nu$ and omit the dependencies on the choice
of fit set and the optimization criterion used to determine the fit
coefficients $c_{\mu\nu\alpha}$. The exact representation of
$\hat{\mathcal{K}}$ in $X_{p}$ can then be written as

$\displaystyle\mathcal{K}_{\mu\nu\kappa\lambda}=\sum_{\alpha\beta}c_{\mu\nu\alpha}\mathcal{K}_{\alpha\beta}c_{\kappa\lambda\beta}+$
$\displaystyle\sum_{\alpha}c_{\mu\nu\alpha}\int
d\bm{r}f_{\alpha}(\bm{r})\hat{\mathcal{K}}e_{\kappa\lambda}(\bm{r})+$ (6)
$\displaystyle\sum_{\beta}c_{\kappa\lambda\beta}\int
d\bm{r}e_{\mu\nu}(\bm{r})\hat{\mathcal{K}}f_{\beta}(\bm{r})+\int
d\bm{r}e_{\mu\nu}(\bm{r})\hat{\mathcal{K}}e_{\kappa\lambda}(\bm{r}).$

Using the scalar product notation (3) we also have

$\displaystyle(f_{\alpha}|\chi_{\mu}\chi_{\nu})$
$\displaystyle=\sum_{\beta}\mathcal{K}_{\alpha\beta}c_{\mu\nu\beta}+(f_{\alpha}|e_{\mu\nu})$
(7)

Keeping in mind that the function $e_{\mu\nu}$ depends implicitly on the fit
space and on the kernel used to define the scalar product, we can define this
function to lie entirely in $X_{e}=X_{p}\setminus X_{f}$ so that we have:

$\displaystyle(f_{\alpha}|e_{\mu\nu})=0.$ (8)

Note that integrals over $f_{\alpha}$ and $e_{\mu\nu}$ with other kernels are
in general non-zero. Then, assuming $\mathcal{K}_{\alpha\beta}$ invertible and
considering the symmetry of scalar products,we may write

$c_{\mu\nu\alpha}=\sum_{\beta}(\chi_{\mu}\chi_{\nu}|f_{\beta})(\mathcal{K}^{-1})_{\beta\alpha}.$
(9)

This choice of fit coefficients in (9) can also be viewed as minimizing the
Lagrangian

$\displaystyle\mathcal{L}_{\mu\nu}$ $\displaystyle=(e_{\mu\nu}|e_{\mu\nu})$
(10)
$\displaystyle=({\chi_{\mu}\chi_{\nu}}|\chi_{\mu}\chi_{\nu})-2\sum_{\beta}c_{\mu\nu\beta}(\chi_{\mu}\chi_{\nu}|f_{\beta})+\sum_{\alpha\beta}c_{\mu\nu\alpha}c_{\mu\nu\beta}(f_{\alpha}|f_{\beta})$

for every $e_{\mu\nu}$. Minimizing $\mathcal{L}_{\mu\nu}$ guarantees a
reasonable precision also for off-diagonal terms in the residuals matrix since
due to the Cauchy-Schwartz inequality we have

$\displaystyle(e_{\mu\nu}|e_{\kappa\lambda})\leq\sqrt{\mathcal{L}_{\mu\nu}\mathcal{L}_{\kappa\lambda}}\;.$
(11)

With the fit coefficient definition (8), the cross error terms in (6) vanish.
This property arises naturally because the metric used in fitting is defined
through the same kernel we aim to fit. Using $\mathcal{K}=v_{c}$ and functions
of the form (2), the theoretical framework just presented is known as (global)
robust Density Fitting (DF).61, 62, 63, 116, 117, 64 It has the advantage of
reducing the storage complexity of the matrix elements of the kernel and the
amount of integrals to be evaluated from $\mathcal{O}\left(N^{4}\right)$ to
$\mathcal{O}\left(N^{3}\right)$. DF is not an approximation if the expansion
is complete, and in this case a compression would only be achieved for exact
linear dependencies in $X_{p}$. In practice the compression is obtained at the
price of an approximation since for reasons of computational efficiency
$F_{f}$ does not span the complete space of products of primary basis
functions. Considering what is left out, we note that a product set defined as
$\chi_{\mu}\chi_{\nu}$ is strongly non-orthogonal. Orthogonalization of such a
basis to span as much as possible of the full Hilbert space would result in
linear combinations of $\chi_{\mu}\chi_{\nu}$ with large coefficients. Given
the finite precision of computer operations, the calculation of matrix
representations of these parts of $\mathcal{H^{K}}$ is likely to be
numerically unstable. In addition, we can expect such combinations of pair
functions to be of minor physical relevance for a quantum chemical
calculation. For this reason it is also numerically favourable to work with a
fit set that is better behaved in terms of orthogonality than an
orthogonalized product set.

DF reduces the asymptotic scaling of the evaluation of RPA and direct MP2
correlation energies from $\mathcal{O}\left(N^{6}\right)$ and
$\mathcal{O}\left(N^{5}\right)$, respectively to
$\mathcal{O}\left(N^{4}\right)$.60 However, the asymptotic scaling of methods
involving exchange terms is not automatically reduced. For instance, using
(6), the exchange contribution to the Fock matrix can be expressed as118

$F_{\mu\nu}=\sum_{\kappa\lambda}P_{\kappa\lambda}\sum_{\alpha\beta}c_{\mu\kappa\alpha}v_{\alpha\beta}c_{\nu\lambda\beta}\;,$
(12)

which is evaluated with the same asymptotic scaling of
$\mathcal{O}\left(N^{4}\right)$ as the variant without density fitting when no
further approximations are made. This is also generally true for post-HF
methods where exchange terms profit less from the compression of 4-index
tensors.72, 83 Therefore, it is advantageous to eliminate the exchange terms
entirely and introduce empirical scaling factors, as for instance in the SOS-
MP296 or SOS-CC2 methods.119

One can however also greatly improve upon the efficiency of global DF by
constraining eq. (5) in such a way that the number of non-zero-elements in $c$
only grows linearly with system size. For instance, instead of using the
Coulomb kernel directly61, 63 one can avoid to define the scalar product on
the operator $\hat{\mathcal{K}}$ and define another metric with more suitable
properties. This could for instance be a local kernel, like the overlap
kernel63 (also known as RI-SVS) or an attenuated Coulomb kernel.64, 120 These
kernels have been used successfully to lower the complexity of for instance
$GW$121, 87, RPA76, 77, MP284, 85 and CC2122 calculations. The price to be
paid is the loss of robustness, equation (8) is then not fulfilled so that the
cross terms in (6) become non-zero. An alternative approach which introduces
the desired sparsity in the fit-coefficient tensor more directly is PADF. In
PADF the density fit is restricted to pairwise sums only and subsequently
distance cut-offs are introduced. Using Latin uppercase superscripts to denote
the atomic centers of functions, the PADF expansion of products of basis
functions is

$\displaystyle\chi^{A}_{\mu}(\bm{r})\chi^{B}_{\nu}(\bm{r})=\sum_{\alpha\in
A}c_{\nu\mu\alpha}f_{\alpha}(\bm{r})+\sum_{\beta\in
B}c_{\mu\nu\beta}f_{\beta}(\bm{r})+e_{\mu\nu}(\bm{r})\;,$ (13)

replacing the simpler expansion (5). The notation $\alpha\in A$ indicates that
the summation is restricted to fit functions centered on atom $A$.

###  2.2 Fit Set Generation

####  2.2.1 Fit sets from products of basis functions

It is easily understood that the choice of $F$ is of key importance in a PADF
code. Ideally, the fit set should be generated on-the-fly, tailored to the
primary basis at hand and the precision of the expansion (13) should be
adjustable in a systematic way using only a single parameter. Many algorithms
for this task have been developed for global DF.123, 124, 125, 126 An
alternative way to generate fit sets on-the-fly is Cholesky decomposition,127,
128, 129 but this approach is not straightforwardly generalized to codes which
can not evaluate 3-center integrals involving the Coulomb potential
analytically. We are only aware of two algorithms which have been developed
specifically for PADF.94, 130 We here adopt the one by Ihring and coworkers94
which we recapitulate for completeness.

From all unique combinations of AOs centered on atom $A$ (denoted by $X_{A}$)
we build an atom-specific trial fit set $\tilde{F}_{A}$ of functions of form
(2),

$\displaystyle\tilde{F}_{A}=\\{\tilde{f}_{\mu\nu}(\textbf{r})\text{ such that
}$ $\displaystyle
R_{\mu\nu}(|\textbf{r}_{A}|)=R_{\mu}(|\textbf{r}_{A}|)R_{\nu}(|\textbf{r}_{A}|)\text{
and }$ (14) $\displaystyle|l_{\mu}-l_{\nu}|<l_{\mu\nu}<|l_{\mu}+l_{\nu}|\text{
}\forall\text{ }\chi_{\mu}(\textbf{r})\text{ and }\chi_{\nu}(\textbf{r})\in
X_{A}\\}$

We then regroup $F_{A}$ in subsets with same angular momentum

$\displaystyle\tilde{F_{A}}=\bigcup_{l,m}\tilde{F}_{A,lm}=\bigcup_{l,m}{Y}_{l}^{m}(\theta,\phi)\times\tilde{R}_{A,lm}$
(15)

where $\tilde{R}_{A,lm}$, the set of radial components centered on the same
atom $A$, is multiplied element-wise to the same spherical harmonic. We then
compute the eigenvectors of the matrix
$(\tilde{R}_{\alpha}|\tilde{R}_{\beta})$ for $\tilde{R}_{\alpha}\text{ and
}\tilde{R}_{\beta}\in\tilde{R}_{A,lm}$ and we keep only the ones with
eigenvalue greater than a specific threshold $\epsilon_{\text{fit}}$. This
threshold can be seen as a parameter tuning the fit quality. Setting it to
zero does not imply exact fitting, as it only solves the one-center part
exactly. In addition, as mentioned above, choosing a too small parameter will
likely introduce numerical instabilities. The set of remaining eigenvectors is
called $R_{A,lm}$. Our final fit set is then

$\displaystyle
F=\bigcup_{Alm}F_{A,lm}=\bigcup_{Alm}Y_{l}^{m}(\theta,\phi)\times R_{A,lm}$
(16)

Alg. 1 shows a pseudocode of the algorithm. In alg. 1,
$R_{\mu}(r)\times\\{Y^{m}_{l_{u}}(\theta,\phi)\\}$ denotes basis functions
with same radial part and different spherical harmonic.

Algorithm 1 Fit set generation algorithm

1:for every atom $A$ do

2: for every $R_{\mu}(r)\\{Y^{m}_{l_{u}}(\theta,\phi)\\}$ in $X_{A}$ do

3: for every $R_{\nu}(r)\\{Y^{m}_{l_{v}}(\theta,\phi)\\}$ in $X_{A}$ do

4: $\tilde{F}_{A}\texttt{ push back
}R_{\mu}(r)R_{\nu}(r)\\{Y^{-l}_{l=|l_{u}-l_{v}|},...,Y^{+l}_{l=|l_{u}+l_{v}|}\\}$

5: for every $\tilde{F}_{A,lm}$ in $\tilde{F}_{A}$ do

6: for every couple $\tilde{f}_{a},\tilde{f}_{b}$ in $\tilde{F}_{A,lm}$ do

7: $V_{ab}=(\tilde{f}_{a}(\textbf{r})|\tilde{f}_{b}(\textbf{r}))$

8: Find Spectrum of $V\Rightarrow
V_{ab}=v_{a\alpha}\varepsilon_{\alpha\beta}v^{\dagger}_{\beta b}$ (eigenvalues
in increasing order)

9: $P_{\alpha b}=\sqrt{\varepsilon_{\alpha\beta}}v^{\dagger}_{\beta b}$ ( $m$
s.t. $\varepsilon_{\alpha\alpha}>\epsilon_{a}$ $\Rightarrow\alpha=1..m$ )

10: $\tilde{F}_{A,lm}=P\cdot\tilde{F}_{A,lm}$ (sets assumed to be matrices,
same ordering used in V)

For small basis sets, the fit set generated through such a procedure does
sometimes not lead to sufficiently accurate results. To overcome this problem,
Ihrig et al.94 artificially enlarged the fit set. This is achieved by adding a
new function to each $X_{A}$. In our implementation, the function is a Slater
type orbital with angular momentum $l^{A}_{max}+1$ where $l^{A}_{max}$ is the
maximum angular momentum present in $X_{A}$, and the arbitrarily chosen
exponent $\alpha$ equal to its angular momentum.

$\displaystyle\chi(\textbf{r})=r^{\alpha}\exp(-\alpha
r)Y_{\alpha}(\theta,\phi)$ $\displaystyle\alpha=l^{A}_{max}+1$ (17)

In the following, we will refer to this procedure as $L$-enhancement ($L$-e).

####  2.2.2 Slater type orbital fit sets

As an alternative to the algorithm just described, we herein also test the use
of hand-optimized STO type fit sets. These come with the disadvantage that
they are not systematically improvable through the adjustment of a single
parameter $\epsilon_{\text{fit}}$. They are however more compact and therefore
more suitable for large-scale applications. In this work we use three
different thresholds which we refer to as _Normal_ , _Good_ and _VeryGood_.
The former two contain STO type functions with angular momentum up to $l=4$,
while the latter one contains functions with exponents up to $l=6$. We have
described these fit sets in ref. 83 to which we refer for more details.

###  2.3 Projection Methods

The improved algorithmic scaling deriving from PADF comes with downsides which
needs to be handled carefully in order to retain sufficient numerical
precision. Therefore, we use two related projection methods which we describe
in the following:

####  2.3.1 Projection Method for the Basis Set

To prevent instabilities due to near-linear dependencies in the AO basis set
itself, this basis set size is often reduced by modifying the Löwdin
orthornomalization131, 132 step. The Löwdin transformation matrix
$\mathbf{S}^{-1/2}$ to an orthonormal set follows from the eigensystem of the
overlap matrix

$\displaystyle\mathbf{S}=$ $\displaystyle\mathbf{U}\mathbf{D}\mathbf{U}^{T}$
(18) $\displaystyle\mathbf{S}^{-1/2}=$
$\displaystyle\mathbf{U}\mathbf{D}^{-1/2}\mathbf{U}^{T}\;.$ (19)

Here we can choose to ignore eigenvectors with eigenvalues below a threshold
$\epsilon_{\text{bas}}$. The simplest way to achieve this is to set the
corresponding eigenvectors in $U$ to zero (thus introducing artificial
states). This is done in the ADF implementation. In the BAND implementation,
we define a smaller orthonormal basis by introducing the regularized Löwdin
transformation

$\tilde{\mathbf{S}}^{-1/2}=\mathbf{S}^{-1/2}\tilde{\mathbf{U}}$ (20)

where $\tilde{\mathbf{U}}$ is the non square (tall) matrix obtained by
removing the eigenvectors columns with eigenvalues smaller than
$\epsilon_{\text{bas}}$ from $\mathbf{U}$. Doing so, fewer orbitals are
obtained and appearance of artificial states is avoided. The elimination of
problematic orbitals that are expressed in the original basis with large
coefficients also helps to prevent problems later on when considering the
product basis and can therefore be beneficial to mitigate errors resulting
from the density fitting.

####  2.3.2 HF projection method

A known issue of PADF is that contributions to the electron repulsion energies
can become negative which can lead to variational collapses of HF SCF
calculations91, 133. As we will argue below, this problem is to a large extent
due to integrals over product functions that are difficult to describe by the
fit set. A way to avoid this problem would be to Cholesky decompose the matrix
(1) and keep only the most important Cholesky vectors, but this is not
practical in calculations with Slater or numerical type orbitals. Instead, the
PADF implementation in AMS uses a simple projector in the original AO space

$\displaystyle\mathbf{T}=\mathbf{T}\mathbf{R}\mathbf{T}^{T}$ (21)

where we use the eigensystem of the AO overlap matrix

$\displaystyle\mathbf{S}=\mathbf{U}\mathbf{D}\mathbf{U}^{T}$ (22)

with $\mathbf{D}$ being a diagonal matrix with the eigenvalues on the
diagonal, and $\mathbf{U}$ having the eigenvectors stored as columns. The
diagonal matrix $\mathbf{R}$ is obtained from $\mathbf{D}$ by taking

$\displaystyle R_{ij}=\delta_{ij}\Theta(D_{ii}-\epsilon_{K})$ (23)

Note that for $\epsilon_{K}=0$ we get $\mathbf{R}=\mathbf{1}$ and hence
$\mathbf{T}=\mathbf{1}$. Heuristically, the action of $\mathbf{T}$ on a vector
or matrix is to remove components of the space parallel to eigenvectors of the
overlap matrix with eigenvalues smaller than the specified threshold
$\epsilon_{K}$. At the SCF stage, the projector is applied both on the left
and on the right of the exact exchange matrix $\mathbf{K}$

$\displaystyle\tilde{\mathbf{K}}=\mathbf{T}^{T}\mathbf{K}\mathbf{T}$ (24)

While the regularized Löwdin orthonormalization removes a subspace for all
energy terms, the HF projector method neglects only the (small) stabilizing
action of the exchange energy, shifting energies upwards. The default value
for $\epsilon_{K}=10^{-3}$ in BAND can therefore be much higher than the one
for the Löwdin orthonormalization projector $\epsilon_{\text{bas}}=10^{-8}$.

We also use the same projector (21) to calculate correlation energies. For
this we redefine the matrix $\mathbf{b}$ which transforms from the AO to the
MO basis as

$b^{\prime}_{i\mu}=\sum_{\mu^{\prime}}b_{i\mu^{\prime}}T_{\mu^{\prime}\mu}\;.$
(25)

Here, the use of the projector is supposed to improve the accuracy of
correlation energies by removing a subspace leading to AO products which can
only be represented poorly by the fit set. The usefulness of the application
of the PM-$K$ to correlation energies can be rationalized as follows: If we
consider an eigenvector
$v_{o}(\textbf{r})=\sum_{\mu}s_{o\mu}\chi_{\mu}(\textbf{r})$ of the overlap
matrix, we can notice that $\sum_{o}s^{\dagger}_{\mu
o}s_{o\nu}=D^{-1}_{\mu\nu}$. From this we understand that the average order of
magnitude of the coefficients is at least of the order of
$s_{o\mu}\sim\frac{1}{\sqrt{D_{\mu\mu}}}$. This shows that the basis set
poorly describes eigenvectors relative to small eigenvalues of $S$, since
large coefficients are needed to expand a small orthogonal component. We then
expect that products involving such linear combinations are the most difficult
to fit. This is mostly because of the diffuse products of basis functions from
distant atoms and from the consequent difficulties in using the PADF
approximation to express such products.113

###  2.4 RPA and SOS-MP2 correlation energies

We now briefly discuss how PADF can be used to speed up the evaluation of SOS-
MP2 and RPA correlation energies, summarizing the more detailed discussions in
ref. 83 and 113. In the basis of fit functions, the RPA correlation energy can
be expressed as

$\displaystyle E^{\text{RPA}}_{c}=$
$\displaystyle\frac{1}{2\pi}\int^{\infty}_{0}d\omega\text{Tr}\left\\{\left[\log\left(\mathbf{1}-\mathbf{Z}(i\omega)\right)\right]+\mathbf{Z}(i\omega)\right\\}$
(26) $\displaystyle=$
$\displaystyle\sum^{N_{\omega}}_{k}\sigma_{k}\text{Tr}\left\\{\left[\log\left(\mathbf{1}-\mathbf{Z}(i\omega_{k})\right)\right]+\mathbf{Z}(i\omega_{k})\right\\}\;,$

which follows directly from the corresponding real-space representation of the
RPA.134 The integration is performed over the imaginary frequency axis for
which either modified Gauss-Legendre grids or, more efficiently,135 minimax
grids of size $N_{\omega}$ can be used. $\omega_{k}$, $\sigma_{k}$ denote
points and corresponding integration weights on the imaginary axis.
$\mathbf{Z}$ is

$Z_{\alpha\beta}(i\omega)=\sum_{\gamma}P^{(0)}_{\alpha\gamma,i\omega}v_{\gamma\beta}\;,$
(27)

and is obtained through the non-interacting polarizability $P^{(0)}$ and the
electron-electron interaction $v$ in the basis of fit functions. Since matrix
logarithms are difficult to calculate, we use that (assuming
$\mathbf{1}-\mathbf{Z}$ can be diagonalized with eigenvalues $\lambda_{j}$)

$\text{Tr}\left[\log(\mathbf{1}-\mathbf{Z})\right]=\sum_{j}\log(\lambda_{j})=\log\left(\prod_{j}\lambda_{j}\right)=\log|\mathbf{1}-\mathbf{Z}|\;,$

and evaluate the determinant $|\mathbf{1}-\mathbf{Z}|$ instead. The imaginary
frequency representation of $P^{(0)}$ is obtained from its discrete imaginary
time representation using nodes
$\left\\{\tau_{k}\right\\}_{j=1,\dots,N_{\tau}}$. The transformation is
achieved by the discrete cosine transform (since $P^{(0)}$ is bosonic)

$P^{(0)}(i\omega_{k})=-i\sum^{N_{\tau}}_{j}\gamma^{(c)}_{kj}\cos(\omega_{k}\tau_{j})P^{(0)}(i\tau_{j})\;,$
(28)

where the $\gamma^{(c)}_{kj}$ are the matrix elements of the kernel of the
discrete cosine transform. The imaginary time representation of $P^{(0)}$ is
given as

$P^{(0)}_{\alpha\beta}(i\tau_{j})=-ic_{\mu\nu\alpha}G^{(0),<}_{\mu\kappa}(-i\tau_{j})G^{(0),>}_{\nu\lambda}(i\tau_{j})c_{\kappa\lambda\beta}\;,$
(29)

with greater and lesser components of the time-ordered Green’s functions being
defined as

$\displaystyle G^{(0),<}_{\mu\nu}(i\tau_{j})=$
$\displaystyle-i\sum_{i}b^{\prime}_{\mu
i}e^{-\epsilon_{i}\tau_{j}}b^{\prime}_{i\nu}$ (30) $\displaystyle
G^{(0),>}_{\mu\nu}(i\tau_{j})=$ $\displaystyle-i\sum_{a}b^{\prime}_{\mu
a}e^{-\epsilon_{a}\tau_{j}}b^{\prime}_{a\nu}\;.$ (31)

Notice, that $b^{\prime}$, as defined in (25), appears in these equations so
that the the PM enters the PRA correlation energies through $G^{\lessgtr}$. As
previously discussed in some detail by us,113 the evaluation of (29) scales
asymptotically as $\mathcal{O}\left(N^{2}\right)$ with system size in the PADF
approximation, since the number of elements in $c$ only scales as
$\mathcal{O}\left(N\right)$ with system size. For detailed working equations
we refer to our previous work.83, 113 In practice, the efficiency of this
approach depends on the possibility to represent the imaginary time and
frequency dependencies of the polarizability as compactly as possible. We
follow Kaltak and Kresse135, 136 and use non-uniformly spaced minimax and
least square grids as described in ref. 137 for imaginary time, and ref. 138
for the imaginary frequency domain.

Alternatively, the imaginary frequency polarizability $P^{(0)}(i\omega_{k})$
can be evaluated directly in the MO basis as for instance described in ref.
108,

$P^{(0)}_{\alpha\beta}(i\omega_{k})=-\sum_{a}\sum_{i}c_{ia\alpha}\frac{1}{\epsilon_{a}-\epsilon_{i}-i\omega_{k}}c_{ia\beta}+c.c.$
(32)

where

$c_{ia\alpha}=\sum_{\mu\nu}b^{\prime}_{\mu i}c_{\mu\nu\alpha}b^{\prime}_{\nu
a}\;.$ (33)

Using the series expansion of the logarithm in (26) one obtains the direct
term of the MP2 correlation energy, $E_{c}^{(2)}$ as its second-order term in
$v_{c}$. $E_{c}^{(2)}$ can be evaluated directly in imaginary time and is
given by

$\displaystyle E^{(2)}_{c}=$
$\displaystyle-\frac{1}{2}\sum^{N_{\tau}}_{k}\alpha_{k}\text{Tr}\left(\sum_{\gamma}Z_{\alpha\gamma,\tau_{k}}Z_{\gamma\beta,\tau_{k}}\right)$
(34) $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum^{N_{\tau}}_{k}\alpha_{k}\sum_{\alpha\beta}Z_{\alpha\beta,\tau_{k}}Z_{\beta\alpha,\tau_{k}}\;,$

where $\left\\{\alpha_{k}\right\\}_{k=1,N_{\tau}}$ denote the integration
weights corresponding to the points
$\left\\{\tau_{k}\right\\}_{k=1,N_{\tau}}$. In the spin-polarized case we need
to calculate

$\displaystyle Z_{\alpha\beta,\tau_{k}}Z_{\beta\alpha,\tau_{k}}=$
$\displaystyle\sum_{\sigma=\alpha,\beta}\sum_{\sigma^{\prime}=\alpha,\beta}Z_{\alpha\beta,\sigma,\tau_{k}}Z_{\alpha\beta,\sigma^{\prime},\tau_{k}}\;.$
(35)

When working in the AO basis, we are only interested in the contribution to
$E_{c}^{(2)}$ from electrons with unpaired spins which is used for instance in
spin-opposite scaled (SOS) MP296 or in DOD-DHs.139 In that case, only the
terms with $\sigma\neq\sigma^{\prime}$ contribute and the resulting
correlation energy expression is scaled by an empirical factor.

While PADF can also be used to lower the time complexity of the exchange
contribution to MP2 from $\mathcal{O}\left(N^{5}\right)$ to
$\mathcal{O}\left(N^{3}\right)$83, the resulting working equations can only be
implemented with a very high prefactor and are therefore not useful in
practice. Instead, the full MP2 correlation energy is evaluated in the MO
basis

$\displaystyle
E_{\text{MP2}}=\sum_{ijab}\frac{(ia|jb)^{2}+\frac{1}{2}[(ia|jb)-(ib|ja)]^{2}}{\varepsilon_{i}+\varepsilon_{j}-\varepsilon_{a}-\varepsilon_{b}}\;.$
(36)

The expression is evaluated as in typical DF-MP2 codes.140 The necessary fit-
coefficients are again transformed to the MO basis using the relation (33).

###  2.5 Summary of the projection method thresholds

To improve the numerical accuracy we apply at various stages of the
calculation related techniques which we loosely call projector methods.
Although the eigenvectors of the overlap matrix form an orthogonal basis set,
some of the functions thus constructed are in an Orwellian sense more equal
than others: those with a small eigenvalue are less valuable, and may even be
numerically harmful. In this spirit the first threshold
$\epsilon_{\text{bas}}$ comes from the regularized Löwdin orthonormalization
transformation (20), and is about removing completely a subspace from the
basis set. Next we have the projector (21) associated with an eigenvalue
threshold $\epsilon_{K}$. We can independently apply this to the exchange
matrix $K$ through the similarity tranformation (24) for Hartree Fock
calculations, and to the orbital coefficients like in equation (25) entering
the MP2 and RPA energy expressions. We do not consider the possibility of
using different values for $\epsilon_{K}$ for Hartree Fock, MP2 and RPA, other
than completely bypassing the projector. Finally, when constructing an
automatic fit set for an atom type, again a regularized Löwdin
orthonormalization is used, based on the eigensystem of the overlap matrix of
the fit functions now in the Coulomb metric, controlled by the eigenvalue
threshold $\epsilon_{\text{fit}}$. This last threshold controls the number of
fit functions, see Alg. 1. The same threshold is also employed for the pseudo
inverse to obtain the PADF fit coefficients.

##  3 Computational Details

All calculations have been performed with modified development versions of the
ADF141 and BAND103 modules of AMS2022.97 All Psi4104 calculations have been
performed using version 1.6.1.

###  3.1 AMS calculations

For all calculations using GTOs we used correlation consistent Dunning basis
sets of double-$\zeta$ (DZ), triple-$\zeta$ (TZ), quadruple-$\zeta$ (QZ) and
5Z quality142 from the basis set exchange library.143 For comparison with ADF,
we used the triple-$\zeta$ plus double polarization (TZ2P) basis set.144

In all BAND calculations we used the PADF scheme for the Hartree-Fock exchange
operator while the Hartree potential fitting procedure is based upon a
partitioning of the density in atomic reservoirs each of which is expanded in
products of radial splines and spherical harmonics.145 This procedure does not
rely in any way on the PADF fit functions. Throughout this work we performed
tests varying the threshold $\epsilon_{\text{fit}}$ to control the size of the
fit set. The thresholds used for particular calculations will be indicated in
the next section. The same holds for the threshold for the canonical
orthogonalization of the primary basis. If not indicated otherwise, we set the
numerical quality to _VeryGood_ , which controls the accuracy of the numerical
integration,146 the quality of the ZLMfit,145 as well as of the threshold
controlling distance effects in HF, MP2 and RPA calculations.83 The same
settings have been used in all ADF calculations. In order to make the basis
functions more compact, in BAND the radial part of the basis functions is
multiplied by a Fermi-Dirac (FD) function by default. We disabled this
behavior in all calculations. In all RPA calculations for the S66 dataset we
calculated the polarizability directly in imaginary frequency using (32) and
modified Gauss–Legendre grids as described in ref. 54 with 50 integration
points.

In all calculations for the L7111 and S30L112 datasets we used the imaginary
time based algorithms for RPA and SOS-MP2 algorithm. In all SOS-MP2
calculations we used 12 imaginary time points which ensures $\mu$Hartree
convergence of correlation energies of organic systems with large HOMO-LUMO
gaps.71, 75 In all RPA calculations113, 59 we used 24 imaginary frequency and
imaginary time points each and and used PBE orbitals as input (RPA@PBE).

We calculated the interaction energies of the dimers in the CIM8 dataset68 at
the RPA@PBE level of theory. If not indicated otherwise we used correlation
consistent Dunning basis sets of DZ and TZ quality. We then extrapolated the
final correlation energies to the complete basis set limit using the
relation,147

$E_{\text{CBS}}=E_{xZ}+\frac{E_{xZ}*x^{3}-E_{(x-1)Z}*(x-1)^{3}}{x^{3}-(x-1)^{3}}\;,$
(37)

where $x=3$ for TZ, $x=4$ for QZ, and so on. We set the numerical quality to
_Good_ and set the threshold controlling distance effects for the RPA
calculation to _Normal_. Also here we used various settings for the quality of
the fit set. For details we refer to the next section. For reasons discussed
in the next section, if not indicated otherwise we set $\epsilon_{K}=10^{-2}$
and $\epsilon_{\text{bas}}=5\cdot 10^{-4}$ for calculations on the CIM8
dataset. Detailed input settings for all calculations can be found in the
supporting information.

###  3.2 Psi4 calculations

We performed Psi4 calculations for the S66 database using Dunning
double-$\zeta$ (DZ), triple-$\zeta$ (TZ), quadruple-$\zeta$ (QZ) to
quintuple-$\zeta$ (5Z) basis set, to perform global DF-MP2 calculations (in
the following simply referred to as DF-MP2). We used default settings for all
calculations, We used the default fit sets for each basis set, i.e. cc-pvxZ-
RI118, 148, 149 for the primary basis cc-pvxZ.

##  4 Results

In this section, we assess the accuracy of the algorithms described herein. We
proceed as follows: In section 4.1 we first illustrate the effect of the
threshold chosen for the regularized Löwdin orthogonaliazation (basis set
reduction, $\epsilon_{\text{bas}}$) as well as for the HF projector method
($\epsilon_{K}$), on the exchange matrix, for a simple molecule. In the
subsequent sections, we compare our results for molecules of increasing size,
starting with the S66 database in section 4.2 and moving on to the L7 and S30L
databases which contain molecules with more than 200 atoms in section 4.4.
Finally, in section 4.5 we showcase the capabilities of our PADF based
algorithms by calculating the interaction energies of 7 large non-covalently
bound complexes in the CIM8 set by Neese and coworkers with up to 910 atoms
(4500 electrons).68

###  4.1 Effect of the HF Projector Method

Figure 1: Comparison fit quality $\epsilon_{\text{fit}}$ with PM threshold
$\epsilon_{K}$, in A figure we plot bonding energy of fluorobenzene, in B the
deviation from its MP2 correlation energy computed through Psi4. Picture B
contains a subset of the grid checked in A.

We start the discussion of our results by illustrating the effect of the HF
projector method on the PADF-MP2 correlation energy of the fluorobenzene
molecule for varying size of the fit set. In the heatmaps in Fig. 1A and Fig.
1B, we show the MP2 bonding energy and MP2 correlation energy of fluorobenzene
for different thresholds for the HF PM $\epsilon_{K}$ and fit quality
$\epsilon_{\text{fit}}$. In particular, in 1A we report the bonding energy and
highlight three interesting zones in the heatmap. On the right side, the PM
threshold $\epsilon_{K}$ is large enough to completely remove the exchange
energy thus reducing to the Hartree limit of HF. In the lower left corner
where the quality of the fit set is poor, the PM threshold does only have a
small effect and the PADF-MP2 bonding energy collapses to unphysically low
values. The rest of the heatmap shows a more stable behavior, but still a
broad range of bonding energies. To check the accuracy of the PADF-MP2
implementation in BAND, in Fig. 1B we show values deviating less than 0.1 %
from the reference DF-MP2 correlation energy (Psi4), and we blur in grey the
ones outside of the range. Thus we observe that increasing the threshold
$\epsilon_{K}$, allows to use smaller fit sets while still maintaining a good
precision in the result.

###  4.2 Deviation BAND-Psi4

Figure 2: Absolute differences in non-covalent interaction energies between
Psi4 and BAND and BAND and ADF for the S66 database using different basis
sets. All values expressed in $\nicefrac{{\text{kcal}}}{{\text{mol}}}$

.

Codes | Basis Set | Max. | MAD | Worst molecule
---|---|---|---|---
$\Delta$MP2 | cc-pVTZ | 0.07 | 0.02 | Pentane-Pentane
cc-pVQZ | 0.06 | 0.02 | Pentane-Pentane
cc-pV5Z | 0.05 | 0.01 | Uracil-Cyclopentane
$\Delta$ (ADF-BAND) | TZ2P | 0.05 | 0.01 | AcOH-Uracil
$\Delta$ RPA | cc-pVTZ | 0.14 | 0.05 | Uracil-Uracil ($\pi-\pi$)
Table 1: Maximum deviations and Mean absolute deviations (MAD) between
different PADF and DF implementations of MP2 and RPA for the non-covalent
interaction energies in the S66 database for different basis sets. All values
expressed in $\nicefrac{{\text{kcal}}}{{\text{mol}}}$.: First three rows:
Deviation of DF-MP2 (Psi4) and PADF-MP2 (BAND), fourth row: Deviations of the
PADF-MP2 implementations in ADF and BAND, last row: Deviations of PADF-RPA
(BAND) to DF-RPA (TURBOMOLE) corrected by the frozen core error (see text for
explanation).

We now compare our PADF-MP2 and PADF-RPA results to DF-MP2 and DF-RPA
reference values for the non-covalent interaction energies of the S66
database. We calculated the DF-MP2 reference values using Psi4. As reference
values for PADF-RPA we use the TURBOMOLE results by Furche and coworkers50 who
used the frozen core approximation. Comparing our DF-MP2 results obtained with
Psi4 to the DF-MP2 results of TURBOMOLE, we found the error introduced by the
frozen core approximation to be between 0.0 and 0.1
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$ for relative energies. Therefore, to
allow for a better comparison of our PADF-RPA to the DF-RPA results, we
corrected the latter ones by the frozen core error for MP2, assuming the
impact of the frozen core approximation to be the same for MP2 and RPA. We
stress that this procedure might not remove the frozen core error completely
and therefore the DF-RPA reference values have certainly higher errors bars
than the DF-MP2 ones

In our calculations the relevant thresholds are the ones for the regularized
Löwdin orthonormalization $\epsilon_{\text{bas}}$, the projector method,
$\epsilon_{K}$ and for the size of the fit set, $\epsilon_{\text{fit}}$. For
the MP2 calculations, we used $\epsilon_{\text{bas}}=10^{-8}$ We set the
threshold for the the projector method (PM) to $\epsilon_{K}=10^{-3}$, and we
used $\epsilon_{\text{fit}}=10^{-12}$ except for some of the $5\zeta$
calculations for which we used $10^{-10}$ instead (see supporting information
for details). The $L$-e, was enabled only for $3\zeta$ calculations. This
corresponds to an (unrealistically sized) fit set which is around 15 times
larger than the primary basis.

The absolute deviations (AD) for all PADF-MP2 interaction energies in the S66
database with respect to DF-MP2 are shown in Fig. 2 for Dunning basis sets of
$3\zeta$ to $5\zeta$ quality. In the same plot, we also show the deviations of
the PADF-MP2 implementations in ADF and BAND using STOs. We see that in all
cases, the deviations are smaller than 0.1
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$, irrespective of the basis set. The
small deviations between BAND and Psi4 should primarily be due to errors
introduced by the PADF approximation. Given that the two implementations
differ also in other aspects, we can, however, not exclude the possibility
that differences in other technical parameters might play a role as well, for
instance differences in the definitions of the numerical integration grids.
Also incompleteness of the fit sets used in Psi4 might play a role. The
deviations of BAND to ADF are mostly due to slightly different integration
grids. In any case, for all practical purposes the agreement between the codes
is excellent. Therefore, we did not investigate the precise origin of these
small discrepancies further. The results of these calculations are summarized
in table 1.

In table 1, we show that both maximum error and MAD of PADF-RPA relative to
DF-RPA are about twice as large as for PADF-MP2. The reason for this might be
that the fit errors are more pronounced due to the presence of the higher
powers of $\mathbf{Z}$ (see (27)). There is also a slightly larger uncertainty
in the reference values due to our approach to subtract to the frozen core
errors in the reference calculations. we also assessed the effect of using a
higher threshold of $\epsilon_{\text{bas}}=10^{-5}$. This did however not
change our results at all, demonstrating the numerical stability of our
method.

Figure 3: Differences in MP2 Interaction energies of the S66 dataset with and
without $l$-e in the fit set. All values are in kcal/mol.

The $l$-e procedure is of key importance for both the GTO and STO triple zeta
basis sets. This is illustrated in Fig. 3 where we show the deviations of the
PADF-MP2 interaction energies with and without $l$-e. For the cc-pVTZ basis
set the $l$-e causes differences of the order of $\sim
1.0\nicefrac{{\text{kcal}}}{{\text{mol}}}$ and of $\sim
0.1\nicefrac{{\text{kcal}}}{{\text{mol}}}$ for TZ2P. The larger deviation in
cc-pVTZ from using $l$-e can be traced back to the fact that an STO is added
to a GTO basis set, thus adding functions with a different radial behaviour to
the fit. For the TZ2P basis set the only difference comes from the higher
angular momentum in the basis. The same comparison of PADF-RPA correlation
energies gave similar results.

###  4.3 Comparison of projector methods

Label |  $\epsilon_{K}$ (HF) |  $\epsilon_{K}$ (MP2) | $\epsilon_{\text{bas}}$
---|---|---|---
None | $\times$ | $\times$ | 1e-8
Bas | $\times$ | $\times$ | 1e-4
HF | 1e-3 | $\times$ | 1e-8
MP2 | $\times$ | 1e-3 | 1e-8
Both | 1e-3 | 1e-3 | 1e-8

Figure 4: Absolute differences in non-covalent MP3 interaction energies
between Psi4 and BAND for the S66 database using the cc-pVQZ basis set and
yarying parameters for the different projector methods. Each line in the
figure correspond to a different combination of the three parameters in the
table. Symbol $\times$ means that the respective method has not been used. All
values are in $\nicefrac{{\text{kcal}}}{{\text{mol}}}$.

In Fig. 4 we compare BAND results using several combinations of the above
mentioned thresholds and show absolute deviation with respect to Psi4 for the
S66 dataset. The names for the different parameters combinations are specified
in the table at the bottom of the same figure. We compare results obtained
without any projector(None), without PM-$K$ but using a tighter threshold for
the regularized Löwdin orthonormalization in the basis set (Bas), using PM-$K$
only at the HF or at the MP2 stage, and finally using PM-$K$ in both cases
(mentioned following the order of the table). The best results are obtained
applying PM-$K$ to MP2 only, followed by applying it to both. This shows that
the usage of PM-$K$ at the HF stage does not contribute much to the overall
accuracy and in fact can even worsen it. This supports our hypothesis from
section 2 since correlation energies are dependent on virtual orbitals which
have more nodes and are more diffuse than the occupied ones. Hence, they are
more difficult to represent in terms of an atomic centered basis and
consequently their products more difficult to express in terms of fit
functions. We emphasize however that we use here a very large fit set. When
fewer fit functions are used, PM-$K$ also needs to be used at the SCF stage to
prevent variational collapse. The calculations using a larger value of the
Löwdin orthonormalization threshold lead to improvements with respect to the
None and HF settings since unstable components of the basis are projected out
from all terms of the Fock matrix. This is decisive for the accuracy of MP2
and RPA correlation energies as we have already seen in the comparisons above.

###  4.4 Interaction energies for L7 and S30L

After having demonstrated the excellent agreement of our PADF-MP2/RPA results
with DF-MP2/RPA also for large basis sets up to 5Z quality, we now discuss the
accuracy of PADF-MP2 for the molecules in the L7 and S30L databases, for which
Furche and coworkers have recently published DF-MP2 reference values.50 We
focus on the direct contribution to the MP2 correlation energy only and
calculate SOS-MP2 interaction energies using our imaginary time based PADF-
SOS-MP2 implementation.83 To represent the imaginary time dependence we chose
here 12 integration points in the interval $[0,\infty)$ which is sufficient to
achieve $\mu$H precision in absolute correlation energies. The reference
values by Furche and coworkers have been calculated with the frozen core
approximation.50 We have however seen for the S66 database, that the maximum
error of this approximation is of the order of 0.1
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$ for relative energies and we do not
expect this to change for larger molecules. Finally, Furche and coworkers only
calculated correlation energies using Dunning basis sets while they used basis
sets of the Kahlsruhe type to calculate their HF or KS energies.50 Therefore,
we do not compare the non-covalent interaction energies for these sets but
only the correlation contributions to them.

####  4.4.1 L7

$\epsilon_{K}$ | $l_{\text{max}}$ | Non-covalent interaction energy | $\Delta_{\text{TURBOMOLE}}$
---|---|---|---
$1e-3$ | 6 |  |  |  |  |  |  |  |
System | Ref. | 1e-6 | 1e-8 | 1e-10 | 1e-12 | 1e-6 | 1e-8 | 1e-10 | 1e-12
c2c2PD | $-33.72$ | $-34.92$ | $-34.46$ | $-34.20$ | $-34.22$ | $-1.20$ | $-0.74$ | $-0.48$ | $-0.50$
c3a | $-22.07$ | $-22.96$ | $-22.65$ | $-22.37$ | $-22.53$ | $-0.89$ | $-0.58$ | $-0.30$ | $-0.46$
c3gc | $-37.81$ | $-39.37$ | $-38.81$ | $-38.38$ | $-38.85$ | $-1.56$ | $-1.00$ | $-0.57$ | $-1.04$
cbh | $-11.49$ | $-13.19$ | $-12.61$ | $-12.38$ | $-12.26$ | $-1.70$ | $-1.12$ | $-0.89$ | $-0.77$
gcgc | $-18.44$ | $-19.64$ | $-19.14$ | $-18.92$ | $-22.37$ | $-1.20$ | $-0.70$ | $-0.48$ | $-3.93$
ggg | $-7.96$ | $-8.32$ | $-8.21$ | $-8.13$ | $-8.59$ | $-0.36$ | $-0.25$ | $-0.17$ | $-0.63$
phe | $-5.15$ | $-5.87$ | $-5.63$ | $-5.51$ | – | $-0.72$ | $-0.48$ | $-0.36$ | –
$\epsilon_{K}$ | $l_{\text{max}}$ | Non-covalent interaction energy | $\Delta_{\text{TURBOMOLE}}$
$1e-3$ | 8 |  |  |  |  |  |  |  |
System | Ref. | 1e-6 | 1e-8 | 1e-10 | 1e-12 | 1e-6 | 1e-8 | 1e-10 | 1e-12
c2c2PD | $-33.72$ | $-34.82$ | $-34.23$ | $-33.75$ |  | $-1.10$ | $-0.51$ | $-0.03$ |
c3a | $-22.07$ | $-22.90$ | $-22.50$ | $-28.58$ |  | $-0.83$ | $-0.43$ | $-6.51$ |
c3gc | $-37.81$ | $-39.26$ | $-38.45$ | $-38.06$ |  | $-1.45$ | $-0.64$ | $-0.25$ |
cbh | $-11.49$ | $-12.96$ | $-12.25$ | $-11.94$ |  | $-1.47$ | $-0.76$ | $-0.45$ |
gcgc | $-18.44$ | $-19.53$ | $-18.97$ | $-29.75$ |  | $-1.09$ | $-0.53$ | $-11.31$ |
ggg | $-7.96$ | $-8.29$ | $-8.15$ | $-8.37$ |  | $-0.33$ | $-0.19$ | $-0.41$ |
phe | $-5.15$ | $-5.82$ | $-5.53$ | $-6.02$ |  | $-0.67$ | $-0.38$ | $-0.87$ |
$\epsilon_{K}$ | $l_{\text{max}}$ | Non-covalent interaction energy | $\Delta_{\text{TURBOMOLE}}$
$5e-3$ | 6 |  |  |  |  |  |  |  |
System | Ref. | 1e-6 | 1e-8 | 1e-10 | 1e-12 | 1e-6 | 1e-8 | 1e-10 | 1e-12
c2c2PD | $-33.72$ | $-34.02$ | $-33.80$ | $-33.69$ | $-33.69$ | $-0.30$ | $-0.08$ | $0.03$ | $0.03$
c3a | $-22.07$ | $-22.15$ | $-22.04$ | $-21.96$ | $-21.97$ | $-0.08$ | $0.03$ | $0.11$ | $0.10$
c3gc | $-37.81$ | $-37.95$ | $-37.76$ | $-37.63$ | $-37.67$ | $-0.14$ | $0.05$ | $0.18$ | $0.14$
cbh | $-11.49$ | $-12.02$ | $-11.85$ | $-11.78$ | $-11.71$ | $-0.53$ | $-0.36$ | $-0.29$ | $-0.22$
gcgc | $-18.44$ | $-18.46$ | $-18.35$ | $-18.29$ | $-18.60$ | $-0.02$ | $0.09$ | $0.15$ | $-0.16$
ggg | $-7.96$ | $-7.95$ | $-7.92$ | $-7.91$ | $-7.93$ | $0.01$ | $0.04$ | $0.05$ | $0.03$
phe | $-5.15$ | $-5.31$ | $-5.22$ | $-5.19$ | $-8.18$ | $-0.16$ | $-0.07$ | $-0.04$ | $-3.03$
Table 2: Comparison of SOS-MP2 contributions to the non-covalent interaction
energies in the L7 database to the ones from Furche and coworkers50 for
different numerical settings using the cc-pvTZ basis set. All values are in
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$. The fit sets have been generated
from the basis products using the threshold $\epsilon_{\text{fit}}$

For the fit sets generated from the products of basis functions, the results
of this comparison for the L7 database can be found in table 2. We first focus
on the upper table. The results here have been obtained using a value of
$\epsilon_{K}=10^{-3}$ as threshold for PM-$K$ (the same value as for S66) and
without the $l$-e method. The results for the different values of
$\epsilon_{k}$ controlling the size of the fit set ranging from $10^{-6}$ to
$10^{-12}$ in table 2 show a slow convergence of the relative SOS-MP2
correlation energies to the DF-SOS-MP2 reference values. However, even with
the already very large fit set corresponding to $\epsilon_{k}=10^{-10}$ (the
number of fit functions is around 10 times larger than the number of primary
basis functions), the maximum deviation is still 0.89
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$ and only for the smallest of the
systems (GGG) in L7, the deviation to the TURBOMOLE results reaches an
acceptable value of 0.17 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$. Moreover,
the results for $\epsilon_{k}=10^{-12}$ start to worsen, which can be related
to the occurrence of linear dependencies in the fit set. For the Phe system in
L7 we even obtained a completely unreasonable value for one of the correlation
energies.

Turning to the second table in table 2, we find that the $l$-e reduces the
error with respect to DF-SOS-MP2 compared to the value obtained with the same
thresholds. However, we already observe numerical instabilities for fit sets
corresponding to $\epsilon_{k}=10^{-10}$. Notice however, that the total size
of this fit set is already larger than the one corresponding to
$\epsilon_{k}=10^{-12}$ without the $l$-e. Therefore, also the $l$-e method
does not resolve the issues of PADF for the molecules in L7.

We now turn to the third table in table 2 which shows results obtained without
the $L$-e method but with a larger threshold of the HF projector method,
$\epsilon_{K}=5\cdot 10^{-3}$, instead of $\epsilon_{K}=10^{-3}$ which has
been used in the two previous tables. This changes the PADF-SOS-MP2
correlation energies drastically and bring them into much better agreement
with the TURBOMOLE results. Already for a moderate size of the fit set
corresponding to $\epsilon_{k}=10^{-6}$, the maximum deviation is reduced to
0.53 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$ and for a value of
$\epsilon_{k}=10^{-8}$, for 6 out of 7 systems the agreement with the
reference values is better than 0.1 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$.
Notice again, that 0.1 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$ is of the
order of uncertainty in the reference values due to the frozen core
approximation. The CBH complex is the only system for which the deviation to
the reference is still relatively large. Only with $\epsilon_{k}=10^{-12}$,
the deviation to the TURBOMOLE results reduces to an acceptable value of
$0.22\nicefrac{{kcal}}{{mol}}$. However, first this fit set is very large and
therefore not very useful in applications to large molecules, and second, we
can still observe quite large deviations to TURBOMOLE data for other systems,
most drastically for PHE.

$\epsilon_{K}=5\times 10^{-3}$ | SOS-MP2 | RPA
---|---|---
System | $E^{\text{SOS-MP2}}_{\text{Ref.}}$ | normal | good | vg | $E^{\text{RPA}}_{\text{Ref.}}$ | normal | good | vg
c2c2PD | $-33.72$ | $-2.09$ | $-0.41$ | $0.09$ | $-30.29$ | $-2.17$ | $-0.32$ | $0.30$
c3a | $-22.07$ | $-1.32$ | $-0.20$ | $0.14$ | $-21.21$ | $-1.43$ | $-0.12$ | $0.32$
c3gc | $-37.81$ | $-2.16$ | $-0.30$ | $0.27$ | $-39.89$ | $-0.19$ | $2.02$ | $2.76$
cbh | $-11.49$ | $-0.83$ | $-0.10$ | $-0.03$ | $-16.84$ | $-0.94$ | $-0.03$ | $0.10$
gcgc | $-18.44$ | $-0.72$ | $-0.05$ | $0.21$ | $-22.33$ | $-0.89$ | $0.04$ | $0.37$
ggg | $-7.96$ | $-0.28$ | $-0.03$ | $0.07$ | $-9.00$ | $-0.32$ | $0.02$ | $0.15$
phe | $-5.15$ | $-0.51$ | $-0.04$ | $0.02$ | $-8.53$ | $-0.60$ | $-0.01$ | $0.10$
$\epsilon_{K}=10^{-2}$ | SOS-MP2 | RPA
System | $E^{\text{SOS-MP2}}_{\text{Ref.}}$ | normal | good | vg | $E^{\text{RPA}}_{\text{Ref.}}$ | normal | good | vg
c2c2PD | $-33.72$ | $-1.31$ | $-0.01$ | $0.09$ | $-30.29$ | $-1.24$ | $0.22$ | $0.71$
c3a | $-22.07$ | $-0.84$ | $0.06$ | $0.14$ | $-21.21$ | $-0.76$ | $0.28$ | $0.62$
c3gc | $-37.81$ | $-1.37$ | $0.11$ | $0.27$ | $-39.89$ | $0.90$ | $2.68$ | $3.25$
cbh | $-11.49$ | $-0.73$ | $-0.07$ | $-0.03$ | $-16.84$ | $-0.73$ | $0.09$ | $0.17$
gcgc | $-18.44$ | $-0.57$ | $0.07$ | $0.21$ | $-22.33$ | $-0.66$ | $0.28$ | $0.55$
ggg | $-7.96$ | $-0.17$ | $0.02$ | $0.07$ | $-9.00$ | $-0.15$ | $0.12$ | $0.22$
phe | $-5.15$ | $-0.43$ | $-0.05$ | $0.02$ | $-8.53$ | $-0.43$ | $0.04$ | $0.11$
Table 3: Deviations $\Delta_{\text{TURBOMOLE}}$ of SOS-MP2 (left) and RPA
(right) contributions to the non-covalent interaction energies in the L7
database to the ones from Furche and coworkers50 for different numerical
settings using the cc-pvTZ basis set. All values are in
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$. Standard STO type fit sets of
varying size, randing from _Basic_ to _VeryGood_ (vg) quality have been used.

##### STO type fit sets

Out of all systems in L7, CBH contains the largest number of atoms, but it is
not the system with most electrons. Therefore, it can be considered as the
most spatially extended system. As we have already discussed, such systems are
expected to be most problematic for PADF based methods since many products of
diffuse functions will occur. Since GTOs decay much faster than STOs from
their atomic centres, it is natural to ask whether GTOs are the best choice to
fit such products of diffuse atomic orbitals. Therefore, we investigate the
accuracy which can be achieved by fitting the products of GTOs with STO type
functions in table 3. In the second table table of table 3 we show some
results for $\epsilon_{K}=10^{-2}$. Additionally, we have reduced the
integration quality to _Good_ and the quality of the threshold controlling
distance effects in the SOS-MP2 method to _Normal_. Both settings together
drastically speed up the PADF-SOS-MP2 calculations.

Using the _Normal_ fit set, relatively large errors are obtained. The results
using the _Good_ and _VeryGood_ fit sets are in relatively good agreement with
TURBOMOLE. Especially for the CBH complex, the deviation is much smaller than
for the GTO type fit sets. Increasing $\epsilon_{K}$ to $10^{-1}$ improves
agreement with the TURBOMOLE results further and using the _Good_ fit set
results in perfect agreement with the reference values (given their
uncertainties due to the frozen core approximation) are obtained.

Independently of the value of $\epsilon_{K}$, the errors of the relative RPA
correlation energies are of the same order of magnitude as for SOS-MP2.
However, as already observed for the S66 benchmark, the errors tend to be
slightly larger. We also observe that relative energies tend to be less
negative, indicating that smaller values of $\epsilon_{K}$ than for PADF-SOS-
MP2 are beneficial for PADF-RPA results. Also here, using the _Normal_ fit
set, errors hardly exceed 1 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$. The
system c3gc (A GC base pair absorbed on a Circumcoronene molecule) is
particularly problematic with errors of the order of 3
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$ for the _VeryGood_ fit set. We show
in the supporting information that lowering of $\epsilon_{K}$ leads to much
better agreement with experiment. This shows that the optimal settings for the
projector methods are system specific and further research is needed to
understand these patterns.

##### PADF-SOS-MP2 results using the CC-pVQZ basis set

|  | Non-covalent interaction energy | $\Delta_{\text{TURBOMOLE}}$
---|---|---|---
System | Ref. | 1e-6 | 1e-8 | vg | 1e-6 | 1e-8 | vg
c2c2PD | $-35.47$ | $-35.63$ |  | $-36.66$ | $-0.16$ |  | $-1.19$
c3a | $-23.43$ | $-23.54$ |  | $-24.15$ | $-0.11$ |  | $-0.72$
c3gc | $-40.33$ | $-40.58$ |  | $-41.50$ | $-0.25$ |  | $-1.17$
cbh | $-12.27$ | $-12.64$ | $-12.43$ | $-12.74$ | $-0.37$ | $-0.16$ | $-0.47$
gcgc | $-20.03$ | $-20.22$ |  | $-20.53$ | $-0.19$ |  | $-0.50$
ggg | $-8.48$ | $-8.53$ |  | $-8.67$ | $-0.05$ |  | $-0.19$
phe | $-6.40$ | $-6.62$ | $-6.53$ | $-6.74$ | $-0.22$ | $-0.13$ | $-0.34$
Table 4: Comparison of SOS-MP2 contributions to the non-covalent interaction
energies in the L7 database to the ones from Furche and coworkers50 for
different numerical settings using the CC-pvQZ basis set. All values are in
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$. All values are obtained using
$\epsilon_{K}=5\cdot 10^{-3}$.

Lastly, we examine the quality of the PADF approximation at the QZ level. As
for the S66 database, the results shown in table 4 demonstrate that the
quality of the interaction energies is not deteriorated compared to the TZ
level. Already with the rather moderate threshold of $\epsilon_{k}=10^{-6}$,
the maximum deviation is -0.37 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$ for
CBH. This value is reduced to -0.16 $\nicefrac{{\text{kcal}}}{{\text{mol}}}$
when the threshold is decreased to $\epsilon_{k}=10^{-8}$. Only using the
_VeryGood_ Slater type fit set leads to worse results than at the TZ
level.This is due to the fact that the _VeryGood_ fit set only contains
functions with $l\leq 6$, with the number of $l=6$ functions being rather
small, while the product basis contains functions with $l\leq 8$, with a large
number of functions with $l\geq 6$.

####  4.4.2 S30L

Figure 5: Deviations of PADF-SOS-MP2 and PADF-RPA with different STO type fit
sets with respect to the TURBOMOLE reference values using global DF. All
calculations have been performed using the cc-pvTZ basis set and the settings
from table 3. All deviations are in kcal/mol.

To further assess the quality of PADF-SOS-MP2, we also calculate the SOS-MP2
and RPA contributions to the interaction energies in the S30L set, which
contains 30 non-covalently bound complexes with up to 205 atoms, out of which
8 are charged.We have omitted the systems 15 and 16 which contain Iodine atoms
and for which a comparison to TURBOMOLE would be difficult due to the use of
pseudopotentials in ref. 50. We also do not show results for 3 of the charged
systems due to convergence difficulties. The deviations to the respective
TURBOMOLE results are shown in figure 5 and confirm the observations for the
L7 database in table 3. When the _Good_ fit set is used, the SOS-MP2 and RPA
errors never exceed 0.5 and 1.0 kcal/mol, respectively. These results also
demonstrate that the errors in relative energies do not seem to depend much on
the sizes of the complexes any more after a certain system size is reached.
This is in fact the expected results, since AOs centered on two atoms very far
apart form each other will not overlap and therefore not contribute to the
PADF errors.

###  4.5 RPA Interaction energies of large complexes

After having demonstrated the relatively good accuracy of the PADF
approximation also for the correlation energies of large molecules, we now
calculate the interaction energies in the CIM (Cluster-in-molecule)8 set by
Neese and coworkers68 at the RPA@PBE level of theory. It comprises 8 large
non-covalently bound complexes ranging in size from 200 to 1027 atoms, with
interactions dominated by either $\sigma$-$\sigma$ dispersion or hydrogen
bonding.68

Given the size of these systems, it is clear that these interaction energies
can only be calculated if certain approximations are introduced. Neese and
coworkers used an approach in which they decomposed the complex into smaller
clusters and calculated the correlation energies using a cluster expansion of
the general form

$E^{\text{corr}}=\sum_{I}E^{\text{corr}}_{I}+\sum_{I<J}E^{\text{corr}}_{IJ}\;,$
(38)

where $I$ denotes a subset of localized occupied and virtual molecular
orbitals.68 They proposed to evaluate the correlation energy by two RI-MP2
calculations using the (aug-)cc-pvDZ and the (aug-)cc-pVTZ basis sets and by a
DLPNO-CCSD(T) correction using the smaller basis. They used the augmented
basis sets only for the two smallest systems while for the other systems they
used the non-augmented basis sets.

Basis set incompleteness errors aside, there are major sources of inaccuracies
with this approach, which could potentially lead to large errors. First, both
the CIM and the DLPNO approximations can introduce errors in relative energies
of several $\nicefrac{{\text{kcal}}}{{\text{mol}}}$ for non-covalent
interactions, even when ”tight” truncation thresholds for the DLPNO settings
are used.28, 29 Second, the chosen extrapolation scheme assumes that the CC
basis set error can be faithfully estimated at the MP2 level. This is however
not necessarily the case since the major basis set errors in CC calculations
arise from the direct MP2 contribution,151, 152, 153 a observation termed as
interference effect by Petersson et al.154 Especially since MP2 correlation
energies will be rather inaccurate for the large molecules in CIM8, the
strategy to estimate the basis set incompleteness error at the MP2 level will
be error-prone. Finally, we mention the recently observed disagreement of well
converged CCSD(T) interaction energies with quantum diffusion Monte Carlo
methods for large non-covalently bound complexes.26 For all these reasons, the
values by Neese and coworkers are certainly not of quantitative accuracy.
Despite all this, they are the most accurate reference values which are
available for these large systems and certainly serve as useful frame of
reference for our RPA@PBE calculations.

We have calculated all RPA@PBE interaction energies using (37) using cc-pvDZ
and cc-pVTZ basis sets with and without counterpoise corrections calculations.
We have extrapolated the RPA correlation energies only, while the TZ results
has been used for the remaining components of the interaction energies. We
have verified the accuracy of this strategy by comparison to results using a
(T,Q) extrapolation with the cc-pVTZ and cc-pVQZ basis sets. As shown in table
S1 in the supporting information, the results of both extrapolation schemes
differ by about 6.5 kcal/mol. This also indicates that for systems as large as
the ones in CIM8 for which QZ calculations are out of reach, basis set
incompleteness errors are typically much larger than the errors introduced by
the PADF approximation.

|  |  |  $E_{RPA}$ (D,T) (% cp) |  |  $\Delta E$ (% cp)
---|---|---|---|---|---
System | $N_{\text{atom}}$ | $N_{\text{bas}}$ | 100 % | 0 % | E(ref.) | 100 % | 0 %
1 | $200$ | $5360$ | $-58.22$ | $-67.80$ | $-70.11$ | $11.89$ | $2.31$
2 | $296$ | $6576$ | $-55.76$ | $-62.79$ | $-63.61$ | $7.85$ | $0.82$
4 | $328$ | $9072$ | $-31.01$ | $-34.67$ | $-36.55$ | $5.54$ | $1.88$
3 | $381$ | $10\,806$ | $-14.09$ | $-24.90$ | $-17.83$ | $3.74$ | $-7.07$
5 | $552$ | $12\,080$ | $-34.31$ | $-48.67$ | $-40.13$ | $5.82$ | $-8.54$
6 | $750$ | $17\,316$ | $-58.97$ | $-86.81$ | $-78.80$ | $19.83$ | $-8.01$
7 | $910$ | $21\,932$ | $-336.83$ | $-414.75$ | $-416.08$ | $79.25$ | $1.33$
8 | $1027$ | $22\,778$ | $-25.58$ | $-41.27$ | $-35.70$ | $10.12$ | $-5.57$
MAD: |  |  |  |  |  | $18.00$ | $4.44$
Table 5: Interaction energies for eight large-non-covalently bound complexes.
The number of basis functions refers to the full complex using cc-pVTZ. The
ref. energy in the second last column has been taken from Neese and
coworkers68 and has been calculated at the CIM-DLPNO-CCSD(T)—CIM-RI-MP2(D,T)Z
level of theory.68 All interaction energies are in
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$ and have been extrapolated using (37)
with cc-pVDZ and cc-pVTZ. The counterpoise (cp) correction used is indicated
by a percentage (0% is no cp).

The results of our RPA@PBE calculations are shown in table 5. The counterpoise
corrected results are the most accurate RPA@PBE interaction energies we can
calculate for these large systems. Given the good accuracy of RPA@PBE for
large non-covalently bound complexes50, they might serve as reference values
for more approximate methods. However, they can not be compared directly to
the values by Neese and coworkers since they did not correct for basis set
superposition errors. For this purpose, we also calculated interaction
energies which are not counterpoise corrected.

Overall, reasonable agreement of our RPA@PBE results with the reference values
by Neese and coworkers68 is observed. With a MAD of 4.44
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$, the deviations are of the same order
of magnitude than popular (dispersion corrected) density functionals155, 68,
156 (e.g. $\omega$B97-X-D with a MAD of 5.06
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$ or B3LYP-D4 with a MAD of 4.81
$\nicefrac{{\text{kcal}}}{{\text{mol}}}$) and smaller than the ones for other
_ab initio_ methods like SCS-MP2.68 We emphasize again that even though the
agreement to the CIM-DLPNO-CCSD(T) reference is satisfactory, these non-
counterpoise corrected interaction energies come with large basis set
superposition errors and are therefore likely to be incorrect (see also table
S2 in the supporting information).

##  5 Conclusion

By comparison to DF-MP2 and DF-RPA, we demonstrated the accuracy of PADF-MP2
and PADF-RPA for the S66, L7, and S30L sets of non-covalently bound complexes
ranging from 6 to more than 200 atoms in size. Especially for the small to
medium molecules in the S66 database, PADF comes with negligible loss of
accuracy compared to global density fitting.

The main advantage of PADF over global DF is that is leads to very fast
algorithms for RPA and SOS-MP2. We have shown that the PADF approach is
suitable to calculate the interaction energies of large molecules. In
particular, we calculated the PADF-RPA@PBE interaction energies of eight large
non-covalently bound complexes at the cc-pVTZ level with more than 1000 atoms
and more than 20000 AOs in size on a single compute node.

The choice of fit set is decisive for precise correlation energies with PADF.
We tested two different types of fit sets. In a first variant, the fit set is
generated directly from products of basis functions.94 In a second variant,
fit sets consisting of even tempered series of STOs are used.83 Despite being
much smaller, we found the second kind of fit set to be suitable to express
products of GTOs. This might be due to the slow decay of radial part of the
STOs making them more suitable to fit delocalized products of AOs.

To improve the precision of the PADF approach further, we introduced a
projector which acts directly on the Fock matrix and removes the attractive
component of the excact exchange. most importantly, we also use this method to
project out subspaces of AOs from the orbital coefficients matrix which can
only be represented poorly by the fit set. Especially when smaller fit sets
are used, we showed the PM-$K$ to be of key importance for accurate
interaction energies.

While the PM reduces the correlation energy errors arising from PADF it cannot
completely eliminate them. This is especially true for large molecules, for
which a compromise between accuracy and computational efficiency is required.
For very large systems like the molecules in the CIM8 dataset, it becomes
mandatory to use smaller fit sets which might introduce errors in interaction
energies which can exceed 1 $\nicefrac{{kcal}}{{mol}}$. This is however also
true for many approximations to high-level methods for the calculation of
correlation energies, for instance CC methods based on localized orbitals
using the DLPNO28, 29 or LNO approximations.26 At the moment, the HF projector
method is not used in a system specific way. A computationally efficient way
to do so could be to check the definiteness of the Hartree-exchange matrix at
runtime and to use this information to adjust $\epsilon_{K}$ during the SCF.

In practice, the errors stemming from the PADF approximation will often be of
only minor relevance. Especially MP2 is typically not used by itself but
rather in double hybrids functionals which typically use a fraction of around
30-60 % of MP2 correlation energy,35, 30 scaling the error by the same amount.
Therefore, the already small PADF errors will be negligible for those
functionals. Furthermore, when small fit sets are used, the fit set
incompleteness error always leads to too low interaction energies while basis
set incompleteness errors lead to too high interaction energies. Especially
for large molecules with several hundreds of atoms for which QZ calculations
are not feasible, the fit set incompleteness errors will be much smaller than
basis set incompleteness errors.

On a more general note, it has recently been recognized that methods which
agree well with each other for small and medium molecules might give very
different results for larger systems.26 Also, approximate (dispersion
corrected) GGAs or hybrid functionals which work well for smaller systems are
much more error-prone for large molecules.155, 156 In order to understand the
reasons for this discrepancy, it is mandatory to push the boundaries of first-
principle methods to much larger systems. At the moment, this comes at the
price of numerical errors and further research is needed to develop techniques
to mitigate these errors further.

We thank Erik van Lenthe for fruitful discussions and the idea behind the
projector methods and acknowledge the use of the computing facilities of Vrije
Universiteit Amsterdam and SCM. Edoardo Spadetto acknowledges funding from the
European Union’s Horizon 2020 research and innovation program under grant
agreement No 956813 (2Exciting).

## References

  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , 864–871
  * Kohn and Sham. L. J. 1965 Kohn, W.; Sham. L. J., Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133
  * Levy 1979 Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the ¡i¿v¡/i¿-representability problem. _Proceedings of the National Academy of Sciences_ 1979, _76_ , 6062–6065
  * Hartree 1928 Hartree, D. R. The Wave Mechanics of an Atom with a Non-Coulomb Central Field Part I Theory and Methods. _Math. Proc. Cambridge Philos. Soc._ 1928, _24_ , 89–110
  * Slater 1929 Slater, J. C. The theory of complex spectra. _Phys. Rev._ 1929, _34_ , 1293–1322
  * Fock 1932 Fock, V. Konflgurationsraum und zweite Quantelung. _Zeitschrift für Phys._ 1932, _75_ , 622–647
  * Froese Fischer 1987 Froese Fischer, C. General Hartree-Fock program. _Computer Physics Communications_ 1987, _43_ , 355–365
  * Echenique and Alonso 2007 Echenique, P.; Alonso, J. L. A mathematical and computational review of Hartree–Fock SCF methods in quantum chemistry. _Molecular Physics_ 2007, _105_ , 3057–3098
  * Stöhr et al. 2019 Stöhr, M.; Van Voorhis, T.; Tkatchenko, A. Theory and practice of modeling van der Waals interactions in electronic-structure calculations. _Chem. Soc. Rev._ 2019, _48_ , 4118–4154
  * Burke 2012 Burke, K. Perspective on density functional theory. _The Journal of Chemical Physics_ 2012, _136_ , 150901
  * Reilly and Tkatchenko 2015 Reilly, A. M.; Tkatchenko, A. van der Waals dispersion interactions in molecular materials: beyond pairwise additivity. _Chem. Sci._ 2015, _6_ , 3289–3301
  * Grimme et al. 2016 Grimme, S.; Hansen, A.; Brandenburg, J. G.; Bannwarth, C. Dispersion-Corrected Mean-Field Electronic Structure Methods. _Chem. Rev._ 2016, _116_ , 5105–5154
  * Coester 1958 Coester, F. Bound states of a many-particle system. _Nucl. Phys._ 1958, _7_ , 421–424
  * Coester and Kümmel 1960 Coester, F.; Kümmel, H. Short-range correlations in nuclear wave functions. _Nucl. Phys._ 1960, _17_ , 477–485
  * Čížek 1966 Čížek, J. On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell-Type Expansion Using Quantum-Field Theoretical Methods. _J. Chem. Phys._ 1966, _45_ , 4256–4266
  * Čížek 1969 Čížek, J. On the Use of the Cluster Expansion and the Technique of Diagrams in Calculations of Correlation Effects in Atoms and Molecules. _Adv. Chem. Physics, Vol. XIV_ 1969, _XIV_ , 35–89
  * Paldus et al. 1972 Paldus, J.; Čížek, J.; Shavitt, I. Correlation Problems in Atomic and Molecular Systems. IV. Extended Coupled-Pair Many-Electron Theory and Its Application to the BH3 Molecule. _Phys. Rev. A_ 1972, _5_ , 50–67
  * Raghavachari et al. 1989 Raghavachari, K.; Trucks, G.; John A. Pople,; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. _Chem. Phys. Lett._ 1989, _157_ , 479–483
  * Bartlett and Musiał 2007 Bartlett, R. J.; Musiał, M. Coupled-cluster theory in quantum chemistry. _Rev. Mod. Phys._ 2007, _79_ , 291–352
  * Gyevi-Nagy et al. 2020 Gyevi-Nagy, L.; Kállay, M.; Nagy, P. R. Integral-Direct and Parallel Implementation of the CCSD(T) Method: Algorithmic Developments and Large-Scale Applications. _J. Chem. Theory Comput._ 2020, _16_ , 366–384
  * Pototschnig et al. 2021 Pototschnig, J. V.; Papadopoulos, A.; Lyakh, D. I.; Repisky, M.; Halbert, L.; Severo Pereira Gomes, A.; Jensen, H. J. A.; Visscher, L. Implementation of Relativistic Coupled Cluster Theory for Massively Parallel GPU-Accelerated Computing Architectures. _J. Chem. Theory Comput._ 2021, _17_ , 5509–5529
  * Neese et al. 2009 Neese, F.; Wennmohs, F.; Hansen, A. Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method. _J. Chem. Phys._ 2009, _130_ , 114108
  * Riplinger et al. 2013 Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural triple excitations in local coupled cluster calculations with pair natural orbitals. _J. Chem. Phys._ 2013, _139_ , 134101
  * Nagy and Kállay 2017 Nagy, P. R.; Kállay, M. Optimization of the linear-scaling local natural orbital CCSD(T) method: Redundancy-free triples correction using Laplace transform. _J. Chem. Phys._ 2017, _146_ , 214106
  * Nagy et al. 2018 Nagy, P. R.; Samu, G.; Kállay, M. Optimization of the Linear-Scaling Local Natural Orbital CCSD(T) Method: Improved Algorithm and Benchmark Applications. _J. Chem. Theory Comput._ 2018, _14_ , 4193–4215
  * Al-Hamdani et al. 2021 Al-Hamdani, Y. S.; Nagy, P. R.; Zen, A.; Barton, D.; Kállay, M.; Brandenburg, J. G.; Tkatchenko, A. Interactions between large molecules pose a puzzle for reference quantum mechanical methods. _Nat. Commun._ 2021, _12_ , 1–12
  * Ballesteros et al. 2021 Ballesteros, F.; Dunivan, S.; Lao, K. U. Coupled cluster benchmarks of large noncovalent complexes: The L7 dataset as well as DNA-ellipticine and buckycatcher-fullerene. _J. Chem. Phys._ 2021, _154_ , 154104
  * Brandenburg et al. 2018 Brandenburg, J. G.; Bannwarth, C.; Hansen, A.; Grimme, S. B97-3c: A revised low-cost variant of the B97-D density functional method. _J. Chem. Phys._ 2018, _148_ , 064104
  * Carter-Fenk et al. 2019 Carter-Fenk, K.; Lao, K. U.; Liu, K. Y.; Herbert, J. M. Accurate and Efficient ab Initio Calculations for Supramolecular Complexes: Symmetry-Adapted Perturbation Theory with Many-Body Dispersion. _J. Phys. Chem. Lett._ 2019, _10_ , 2706–2714
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _Journal of Computational Chemistry_ 2006, _27_ , 1787–1799
  * Goerigk and Grimme 2014 Goerigk, L.; Grimme, S. Double-hybrid density functionals. _Wiley Interdiscip. Rev. Comput. Mol. Sci._ 2014, _4_ , 576–600
  * Goerigk et al. 2017 Goerigk, L.; Hansen, A.; Bauer, C.; Ehrlich, S.; Najibi, A.; Grimme, S. A look at the density functional theory zoo with the advanced GMTKN55 database for general main group thermochemistry, kinetics and noncovalent interactions. _Phys. Chem. Chem. Phys._ 2017, _19_ , 32184–32215
  * Mehta et al. 2018 Mehta, N.; Casanova-Páez, M.; Goerigk, L. Semi-empirical or non-empirical double-hybrid density functionals: Which are more robust? _Phys. Chem. Chem. Phys._ 2018, _20_ , 23175–23194
  * Najibi and Goerigk 2018 Najibi, A.; Goerigk, L. The Nonlocal Kernel in van der Waals Density Functionals as an Additive Correction: An Extensive Analysis with Special Emphasis on the B97M-V and $\omega$b97M-V Approaches. _J. Chem. Theory Comput._ 2018, _14_ , 5725–5738
  * Santra et al. 2019 Santra, G.; Sylvetsky, N.; Martin, J. M. Minimally Empirical Double-Hybrid Functionals Trained against the GMTKN55 Database: RevDSD-PBEP86-D4, revDOD-PBE-D4, and DOD-SCAN-D4. _J. Phys. Chem. A_ 2019, _123_ , 5129–5143
  * Møller and Plesset 1934 Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. _Phys. Rev._ 1934, _46_ , 618–622
  * Mezei et al. 2015 Mezei, P. D.; Csonka, G. I.; Ruzsinszky, A.; Kállay, M. Construction and Application of a New Dual-Hybrid Random Phase Approximation. _J. Chem. Theory Comput._ 2015, _11_ , 4615–4626
  * Grimme and Steinmetz 2016 Grimme, S.; Steinmetz, M. A computationally efficient double hybrid density functional based on the random phase approximation. _Phys. Chem. Chem. Phys._ 2016, _18_ , 20926–20937
  * Mezei et al. 2017 Mezei, P. D.; Csonka, G. I.; Ruzsinszky, A.; Kállay, M. Construction of a Spin-Component Scaled Dual-Hybrid Random Phase Approximation. _J. Chem. Theory Comput._ 2017, _13_ , 796–803
  * Macke 1950 Macke, W. Über die Wechselwirkungen im Fermi-Gas. _Zeitschrift fur Naturforsch._ 1950, _5_ , 192–208
  * Langreth and Perdew 1975 Langreth, D. C.; Perdew, J. P. The Exhange-Correlation Energy of a metalic surface. _Solid State Commun._ 1975, _17_ , 1425–1429
  * Heßelmann and Görling 2011 Heßelmann, A.; Görling, A. Random-phase approximation correlation methods for molecules and solids. _Mol. Phys._ 2011, _109_ , 2473–2500
  * Eshuis and Furche 2011 Eshuis, H.; Furche, F. A parameter-free density functional that works for noncovalent interactions. _J. Phys. Chem. Lett._ 2011, _2_ , 983–989
  * Eshuis et al. 2012 Eshuis, H.; Bates, J. E.; Furche, F. Electron correlation methods based on the random phase approximation. _Theor. Chem. Acc._ 2012, _131_ , 1084
  * Ren et al. 2012 Ren, X.; Rinke, P.; Joas, C.; Scheffler, M. Random-phase approximation and its applications in computational chemistry and materials science. _J. Mater. Sci._ 2012, _47_ , 7447–7471
  * Chen et al. 2017 Chen, G. P.; Voora, V. K.; Agee, M. M.; Balasubramani, S. G.; Furche, F. Random-Phase Approximation Methods. _Annu. Rev. Phys. Chem._ 2017, _68_ , 421–445
  * Chedid et al. 2018 Chedid, J.; Ferrara, N. M.; Eshuis, H. Describing transition metal homogeneous catalysis using the random phase approximation. _Theor. Chem. Acc._ 2018, _137:158_ , 1–11
  * Kreppel et al. 2020 Kreppel, A.; Graf, D.; Laqua, H.; Ochsenfeld, C. Range-Separated Density-Functional Theory in Combination with the Random Phase Approximation: An Accuracy Benchmark. _J. Chem. Theory Comput._ 2020, _16_ , 2985–2994
  * Modrzejewski et al. 2020 Modrzejewski, M.; Yourdkhani, S.; Klimeš, J. Random Phase Approximation Applied to Many-Body Noncovalent Systems. _J. Chem. Theory Comput._ 2020, _16_ , 427–442
  * Nguyen et al. 2020 Nguyen, B. D.; Chen, G. P.; Agee, M. M.; Burow, A. M.; Tang, M. P.; Furche, F. Divergence of Many-Body Perturbation Theory for Noncovalent Interactions of Large Molecules. _J. Chem. Theory Comput._ 2020, _16_ , 2258–2273
  * Grüneis et al. 2010 Grüneis, A.; Marsman, M.; Kresse, G. Second-order Møller–Plesset perturbation theory applied to extended systems. II. Structural and energetic properties. _The Journal of Chemical Physics_ 2010, _133_ , 074107
  * Nguyen et al. 2020 Nguyen, B. D.; Chen, G. P.; Agee, M. M.; Burow, A. M.; Tang, M. P.; Furche, F. Divergence of Many-Body Perturbation Theory for Noncovalent Interactions of Large Molecules. _Journal of Chemical Theory and Computation_ 2020, _16_ , 2258–2273, PMID: 32105488
  * Trushin et al. 2021 Trushin, E.; Thierbach, A.; Görling, A. Toward chemical accuracy at low computational cost: Density-functional theory with $\sigma$-functionals for the correlation energy. _J. Chem. Phys._ 2021, _154_ , 014104
  * Fauser et al. 2021 Fauser, S.; Trushin, E.; Neiss, C.; Görling, A. Chemical accuracy with $\sigma$-functionals for the Kohn-Sham correlation energy optimized for different input orbitals and eigenvalues. _J. Chem. Phys._ 2021, _155_ , 134111
  * Erhard et al. 2022 Erhard, J.; Fauser, S.; Trushin, E.; Görling, A. Scaled $\sigma$ -functionals for the Kohn – Sham correlation energy with scaling functions from the homogeneous electron gas. _J. Chem. Phys._ 2022, _157_ , 114105
  * Paier et al. 2012 Paier, J.; Ren, X.; Rinke, P.; Scuseria, G. E.; Grüneis, A.; Kresse, G.; Scheffler, M. Assessment of correlation energies based on the random-phase approximation. _New J. Phys._ 2012, _14_ , 043002
  * Ren et al. 2013 Ren, X.; Rinke, P.; Scuseria, G. E.; Scheffler, M. Renormalized second-order perturbation theory for the electron correlation energy: Concept, implementation, and benchmarks. _Phys. Rev. B - Condens. Matter Mater. Phys._ 2013, _88_ , 035120
  * Hummel et al. 2019 Hummel, F.; Grüneis, A.; Kresse, G.; Ziesche, P. Screened Exchange Corrections to the Random Phase Approximation from Many-Body Perturbation Theory. _J. Chem. Theory Comput._ 2019, _15_ , 3223–3236
  * Förster 2022 Förster, A. Assessment of the second-order statically screened exchange correction to the random phase approximation for correlation energies. _J. Chem. Theory Comput._ 2022, _18_ , 5948–5965
  * Furche 2008 Furche, F. Developing the random phase approximation into a practical post-Kohn-Sham correlation model. _J. Chem. Phys._ 2008, _129_ , 114105
  * Whitten 1973 Whitten, J. L. Coulombic potential energy integrals and approximations. _J. Chem. Phys._ 1973, _58_ , 4496–4501
  * Baerends et al. 1973 Baerends, E. J.; Ellis, D. E.; Ros, P. Self-consistent molecular Hartree—Fock—Slater calculations I. The computational procedure. _Chem. Phys._ 1973, _2_ , 41–51
  * Dunlap et al. 1979 Dunlap, B. I.; Connolly, J. W.; Sabin, J. R. On some approximations in applications of X$\alpha$ theory. _J. Chem. Phys._ 1979, _71_ , 3396–3402
  * Feyereisen et al. 1993 Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Use of approximate integrals in ab initio theory. An application in MP2 energy calculations. _Chem. Phys. Lett._ 1993, _208_ , 359–363
  * Weigend et al. 1998 Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized auxiliary basis sets and demonstration of efficiency. _Chem. Phys. Lett._ 1998, _294_ , 143–152
  * Pinski et al. 2015 Pinski, P.; Riplinger, C.; Valeev, E. F.; Neese, F. Sparse maps - A systematic infrastructure for reduced-scaling electronic structure methods. I. An efficient and simple linear scaling local MP2 method that uses an intermediate basis of pair natural orbitals. _J. Chem. Phys._ 2015, _143_ , 034108
  * Nagy et al. 2016 Nagy, P. R.; Samu, G.; Kállay, M. An Integral-Direct Linear-Scaling Second-Order Møller-Plesset Approach. _J. Chem. Theory Comput._ 2016, _12_ , 4897–4914
  * Ni et al. 2021 Ni, Z.; Guo, Y.; Neese, F.; Li, W.; Li, S. Cluster-in-Molecule Local Correlation Method with an Accurate Distant Pair Correction for Large Systems. _J. Chem. Theory Comput._ 2021, _17_ , 756–766
  * Barca et al. 2021 Barca, G. M.; Galvez Vallejo, J. L.; Poole, D. L.; Alkan, M.; Stocks, R.; Rendell, A. P.; Gordon, M. S. Enabling Large-Scale Correlated Electronic Structure Calculations: Scaling the RI-MP2 Method on Summit. Int. Conf. High Perform. Comput. Networking, Storage Anal. SC. 2021; pp 1–15
  * Zienau et al. 2009 Zienau, J.; Clin, L.; Doser, B.; Ochsenfeld, C. Cholesky-decomposed densities in laplace-based second-order møller-plesset perturbation theory. _J. Chem. Phys._ 2009, _130_ , 204112
  * Doser et al. 2009 Doser, B.; Lambrecht, D. S.; Kussmann, J.; Ochsenfeld, C. Linear-scaling atomic orbital-based second-order Møller-Plesset perturbation theory by rigorous integral screening criteria. _J. Chem. Phys._ 2009, _130_ , 064107
  * Hohenstein et al. 2012 Hohenstein, E. G.; Parrish, R. M.; Martínez, T. J. Tensor hypercontraction density fitting. I. Quartic scaling second- and third-order Møller-Plesset perturbation theory. _J. Chem. Phys._ 2012, _137_ , 1085
  * Hohenstein et al. 2012 Hohenstein, E. G.; Parrish, R. M.; Sherrill, C. D.; Martínez, T. J. Communication: Tensor hypercontraction. III. Least-squares tensor hypercontraction for the determination of correlated wavefunctions. _J. Chem. Phys._ 2012, _137_ , 1101
  * Maurer et al. 2014 Maurer, S. A.; Kussmann, J.; Ochsenfeld, C. Communication: A reduced scaling J-engine based reformulation of SOS-MP2 using graphics processing units. _J. Chem. Phys._ 2014, _141_ , 051106
  * Maurer et al. 2014 Maurer, S. A.; Clin, L.; Ochsenfeld, C. Cholesky-decomposed density MP2 with density fitting: Accurate MP2 and double-hybrid DFT energies for large systems. _J. Chem. Phys._ 2014, _140_ , 224112
  * Wilhelm et al. 2016 Wilhelm, J.; Seewald, P.; Del Ben, M.; Hutter, J. Large-Scale Cubic-Scaling Random Phase Approximation Correlation Energy Calculations Using a Gaussian Basis. _J. Chem. Theory Comput._ 2016, _12_ , 5851–5859
  * Luenser et al. 2017 Luenser, A.; Schurkus, H. F.; Ochsenfeld, C. Vanishing-Overhead Linear-Scaling Random Phase Approximation by Cholesky Decomposition and an Attenuated Coulomb-Metric. _J. Chem. Theory Comput._ 2017, _13_ , 1647–1655
  * Vogler et al. 2017 Vogler, S.; Ludwig, M.; Maurer, M.; Ochsenfeld, C. Low-scaling first-order properties within second-order Møller-Plesset perturbation theory using Cholesky decomposed density matrices. _J. Chem. Phys._ 2017, _147_ , 024101
  * Beuerle et al. 2018 Beuerle, M.; Graf, D.; Schurkus, H. F.; Ochsenfeld, C. Efficient calculation of beyond RPA correlation energies in the dielectric matrix formalism. _J. Chem. Phys._ 2018, _148_ , 204104
  * Graf et al. 2018 Graf, D.; Beuerle, M.; Schurkus, H. F.; Luenser, A.; Savasci, G.; Ochsenfeld, C. Accurate and Efficient Parallel Implementation of an Effective Linear-Scaling Direct Random Phase Approximation Method. _J. Chem. Theory Comput._ 2018, _14_ , 2505–2515
  * Graf et al. 2019 Graf, D.; Beuerle, M.; Ochsenfeld, C. Low-Scaling Self-Consistent Minimization of a Density Matrix Based Random Phase Approximation Method in the Atomic Orbital Space. _J. Chem. Theory Comput._ 2019, _15_ , 4468–4477
  * Duchemin and Blase 2019 Duchemin, I.; Blase, X. Separable resolution-of-the-identity with all-electron Gaussian bases: Application to cubic-scaling RPA. _J. Chem. Phys._ 2019, _150_ , 174120
  * Förster et al. 2020 Förster, A.; Franchini, M.; van Lenthe, E.; Visscher, L. A Quadratic Pair Atomic Resolution of the Identity Based SOS-AO-MP2 Algorithm Using Slater Type Orbitals. _J. Chem. Theory Comput._ 2020, _16_ , 875 – 891
  * Glasbrenner et al. 2020 Glasbrenner, M.; Graf, D.; Ochsenfeld, C. Efficient Reduced-Scaling Second-Order Møller-Plesset Perturbation Theory with Cholesky-Decomposed Densities and an Attenuated Coulomb Metric. _J. Chem. Theory Comput._ 2020, _16_ , 6856–6868
  * Glasbrenner et al. 2021 Glasbrenner, M.; Vogler, S.; Ochsenfeld, C. Efficient low-scaling computation of NMR shieldings at the second-order Møller-Plesset perturbation theory level with Cholesky-decomposed densities and an attenuated Coulomb metric. _J. Chem. Phys._ 2021, _155_ , 224107
  * Drontschenko et al. 2021 Drontschenko, V.; Graf, D.; Laqua, H.; Ochsenfeld, C. Lagrangian-Based Minimal-Overhead Batching Scheme for the Efficient Integral-Direct Evaluation of the RPA Correlation Energy. _J. Chem. Theory Comput._ 2021, _17_ , 5623–5634
  * Wilhelm et al. 2021 Wilhelm, J.; Seewald, P.; Golze, D. Low-scaling GW with benchmark accuracy and application to phosphorene nanosheets. _J. Chem. Theory Comput._ 2021, _17_ , 1662–1677
  * Watson et al. 2003 Watson, M. A.; Handy, N. C.; Cohen, A. J. Density functional calculations, using Slater basis sets, with exact exchange. _J. Chem. Phys._ 2003, _119_ , 6475–6481
  * Manzer et al. 2015 Manzer, S. F.; Epifanovsky, E.; Head-Gordon, M. Efficient Implementation of the Pair Atomic Resolution of the Identity Approximation for Exact Exchange for Hybrid and Range-Separated Density Functionals. _Journal of Chemical Theory and Computation_ 2015, _11_ , 518–527, PMID: 25691831
  * Wirz et al. 2017 Wirz, L. N.; Reine, S. S.; Pedersen, T. B. On Resolution-of-the-Identity Electron Repulsion Integral Approximations and Variational Stability. _J. Chem. Theory Comput._ 2017, _13_ , 4897–4906
  * Merlot et al. 2013 Merlot, P.; Kjaergaard, T.; Helgaker, T.; Lindh, R.; Aquilante, F.; Reine, S.; Pedersen, T. Attractive electron-electron interactions within robust local fitting approximations. _Journal of computational chemistry_ 2013, _34_
  * Hollman et al. 2014 Hollman, D. S.; Schaefer, H. F.; Valeev, E. F. Semi-exact concentric atomic density fitting: Reduced cost and increased accuracy compared to standard density fitting. _J. Chem. Phys._ 2014, _140_ , 064109
  * Hollman et al. 2017 Hollman, D. S.; Schaefer, H. F.; Valeev, E. F. Fast construction of the exchange operator in an atom-centred basis with concentric atomic density fitting. _Mol. Phys._ 2017, _115_ , 2065–2076
  * Ihrig et al. 2015 Ihrig, A. C.; Wieferink, J.; Zhang, I. Y.; Ropo, M.; Ren, X.; Rinke, P.; Scheffler, M.; Blum, V. Accurate localized resolution of identity approach for linear-scaling hybrid density functionals and for many-body perturbation theory. _New J. Phys._ 2015, _17_ , 093020
  * Rebolini et al. 2016 Rebolini, E.; Izsák, R.; Reine, S. S.; Helgaker, T.; Pedersen, T. B. Comparison of Three Efficient Approximate Exact-Exchange Algorithms: The Chain-of-Spheres Algorithm, Pair-Atomic Resolution-of-the-Identity Method, and Auxiliary Density Matrix Method. _J. Chem. Theory Comput._ 2016, _12_ , 3514–3522
  * Jung et al. 2004 Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. Scaled opposite-spin second order moller-plesset correlation energy: An economical electronic structure method. _J. Chem. Phys._ 2004, _121_ , 9793–9802
  * Rüger et al. 2022 Rüger, R.; Franchini, M.; Trnka, T.; Yakovlev, A.; van Lenthe, E.; Philipsen, P.; van Vuren, T.; Klumpers, B.; Soini, T. AMS 2022.1, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, http://www.scm.com. 2022
  * Baerends et al. 2022 Baerends, E.; Ziegler, T.; Atkins, A.; Autschbach, J.; Baseggio, O.; Bashford, D.; Bérces, A.; Bickelhaupt, F.; Bo, C.; Boerrigter, P.; Cavallo, L.; Daul, C.; Chong, D.; Chulhai, D.; Deng, L.; Dickson, R.; Dieterich, J.; Ellis, D.; van Faassen, M.; Fan, L.; Fischer, T.; Förster, A.; Guerra, C. F.; Franchini, M.; Ghysels, A.; Giammona, A.; van Gisbergen, S.; Goez, A.; Götz, A.; Groeneveld, J.; Gritsenko, O.; Grüning, M.; Gusarov, S.; Harris, F.; van den Hoek, P.; Hu, Z.; Jacob, C.; Jacobsen, H.; Jensen, L.; Joubert, L.; Kaminski, J.; van Kessel, G.; König, C.; Kootstra, F.; Kovalenko, A.; Krykunov, M.; van Lenthe, E.; McCormack, D.; Michalak, A.; Mitoraj, M.; Morton, S.; Neugebauer, J.; Nicu, V.; Noodleman, L.; Osinga, V.; Patchkovskii, S.; Pavanello, M.; Peeples, C.; Philipsen, P.; Post, D.; Pye, C.; Ramanantoanina, H.; Ramos, P.; Ravenek, W.; Reimann, M.; Rodríguez, J.; Ros, P.; Rüger, R.; Schipper, P.; Schlüns, D.; van Schoot, H.; Schreckenbach, G.; Seldenthuis, J.; Seth, M.; Snijders, J.; Solà, M.; Stener, M.; Swart, M.; Swerhone, D.; Tognetti, V.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T.; van Wezenbeek, E.; Wiesenekker, G.; Wolff, S.; Woo, T.; Yakovlev, A. ADF2022.1, locally modified development version. 2022
  * Förster and Visscher 2022 Förster, A.; Visscher, L. Quasiparticle Self-Consistent GW-Bethe-Salpeter equation calculations for large chromophoric systems. _J. Chem. Theory Comput._ 2022, _18_ , 6779–6793
  * Wiesenekker et al. 1988 Wiesenekker, G.; te Velde, G.; Baerends, E. J. Analytic quadratic integration over the two-dimensional Brillouin zone. _Journal of Physics C: Solid State Physics_ 1988, _21_ , 4263–4283
  * Wiesenekker and Baerends 1991 Wiesenekker, G.; Baerends, E. J. Quadratic integration over the three-dimensional Brillouin zone. _Journal of Physics: Condensed Matter_ 1991, _3_ , 6721–6742
  * te Velde and Baerends 1991 te Velde, G.; Baerends, E. J. Precise density-functional method for periodic structures. _Phys. Rev. B_ 1991, _44_ , 7888–7903
  * 103 Philipsen, P.; te Velde, G.; Baerends, E.; Berger, J.; de Boeij, P.; Franchini, M.; Groeneveld, J.; Kadantsev, E.; Klooster, R.; Kootstra, F.; Pols, M.; Romaniello, P.; Raupach, M.; Skachkov, D.; Snijders, J.; Verzijl, C.; Gil, J. C.; Thijssen, J. M.; Wiesenekker, G.; Peeples, C. A.; Schreckenbach, G.; Ziegler., T. BAND 2022.1.
  * Turney et al. 2012 Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; Russ, N. J.; Leininger, M. L.; Janssen, C. L.; Seidl, E. T.; Allen, W. D.; Schaefer, H. F.; King, R. A.; Valeev, E. F.; Sherrill, C. D.; Crawford, T. D. Psi4: an open-source ab initio electronic structure program. _WIREs Computational Molecular Science_ 2012, _2_ , 556–565
  * Balasubramani et al. 2020 Balasubramani, S. G.; Chen, G. P.; Coriani, S.; Diedenhofen, M.; Frank, M. S.; Franzke, Y. J.; Furche, F.; Grotjahn, R.; Harding, M. E.; Hättig, C.; Hellweg, A.; Helmich-Paris, B.; Holzer, C.; Huniar, U.; Kaupp, M.; Khah, A. M.; Khani, S. K.; Müller, T.; Mack, F.; Nguyen, B. D.; Parker, S. M.; Perlt, E.; Rappoport, D.; Reiter, K.; Roy, S.; Rückert, M.; Schmitz, G.; Sierka, M.; Tapavicza, E.; Tew, D. P.; van Wüllen, C.; Voora, V. K.; Weigend, F.; Wodynski, A.; Yu, J. M. TURBOMOLE : Modular program suite for ab initio quantum-chemical and condensed-matter simulations. _J. Chem. Phys._ 2020, 184107
  * Blum et al. 2009 Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab initio molecular simulations with numeric atom-centered orbitals. _Comput. Phys. Commun._ 2009, _180_ , 2175–2196
  * Havu et al. 2009 Havu, V.; Blum, V.; Havu, P.; Scheffler, M. Efficient O (N) integration for all-electron electronic structure calculation using numeric basis functions. _J. Comput. Phys._ 2009, _228_ , 8367–8379
  * Ren et al. 2012 Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions. _New J. Phys._ 2012, _14_ , 053020
  * 109 Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. The Fritz Haber Institute ab initio molecular simulations package (FHI-aims), http://www.fhi-berlin.mpg.de/aims, 2009
  * Tew 2018 Tew, D. P. Communication: Quasi-robust local density fitting. _J. Chem. Phys._ 2018, _148_ , 011102
  * Sedlak et al. 2013 Sedlak, R.; Janowski, T.; Pitoňák, M.; Řezáč, J.; Pulay, P.; Hobza, P. Accuracy of quantum chemical methods for large noncovalent complexes. _J. Chem. Theory Comput._ 2013, _9_ , 3364–3374
  * Sure and Grimme 2015 Sure, R.; Grimme, S. Comprehensive Benchmark of Association (Free) Energies of Realistic Host-Guest Complexes. _J. Chem. Theory Comput._ 2015, _11_ , 3785–3801
  * Förster and Visscher 2020 Förster, A.; Visscher, L. Low-Order Scaling G0W0 by Pair Atomic Density Fitting. _J. Chem. Theory Comput._ 2020, _16_ , 7381–7399
  * Řezáč et al. 2011 Řezáč, J.; Riley, K. E.; Hobza, P. S66: A well-balanced database of benchmark interaction energies relevant to biomolecular structures. _J. Chem. Theory Comput._ 2011, _7_ , 2427–2438
  * Loaiza 2017 Loaiza, M. A Short Introduction to Hilbert Space Theory. _Journal of Physics: Conference Series_ 2017, _839_ , 012002
  * Dunlap et al. 1979 Dunlap, B. I.; Connolly, J. W.; Sabin, J. R. On first-row diatomic molecules and local density models. _J. Chem. Phys._ 1979, _71_ , 4993–4999
  * Dunlap and Mei 1983 Dunlap, B. I.; Mei, W. N. Basis set effects on spectroscopic constants for C2 and Si2 and the symmetry dilemma in the X$\alpha$ model. _J. Chem. Phys._ 1983, _78_ , 4997
  * Weigend 2002 Weigend, F. A fully direct RI-HF algorithm: Implementation, optimised auxiliary basis sets, demonstration of accuracy and efficiency. _Phys. Chem. Chem. Phys._ 2002, _4_ , 4285–4291
  * Winter and Hättig 2011 Winter, N. O.; Hättig, C. Scaled opposite-spin CC2 for ground and excited states with fourth order scaling computational costs. _J. Chem. Phys._ 2011, _134_ , 184101
  * Jung et al. 2005 Jung, Y.; Sodt, A.; Gill, P. M. W.; Head-Gordon, M. Auxiliary basis expansions for large-scale electronic structure calculations. _Proc. Natl. Acad. Sci._ 2005, _102_ , 6692–6697
  * Wilhelm et al. 2018 Wilhelm, J.; Golze, D.; Talirz, L.; Hutter, J.; Pignedoli, C. A. Toward GW Calculations on Thousands of Atoms. _J. Phys. Chem. Lett._ 2018, _9_ , 306–312
  * Sacchetta et al. 2022 Sacchetta, F.; Graf, D.; Laqua, H.; Ambroise, M. A.; Kussmann, J.; Dreuw, A.; Ochsenfeld, C. An effective sub-quadratic scaling atomic-orbital reformulation of the scaled opposite-spin RI-CC2 ground-state model using Cholesky-decomposed densities and an attenuated Coulomb metric. _J. Chem. Phys._ 2022, _157_
  * Kállay 2014 Kállay, M. A systematic way for the cost reduction of density fitting methods. _J. Chem. Phys._ 2014, _141_ , 244113
  * Stoychev et al. 2017 Stoychev, G. L.; Auer, A. A.; Neese, F. Automatic Generation of Auxiliary Basis Sets. _J. Chem. Theory Comput._ 2017, _13_ , 554–562
  * Lehtola 2021 Lehtola, S. Straightforward and Accurate Automatic Auxiliary Basis Set Generation for Molecular Calculations with Atomic Orbital Basis Sets. _J. Chem. Theory Comput._ 2021, _17_ , 6886–6900
  * Semidalas and Martin 2022 Semidalas, E.; Martin, J. M. Automatic generation of complementary auxiliary basis sets for explicitly correlated methods. _J. Comput. Chem._ 2022, _43_ , 1690–1700
  * Koch et al. 2003 Koch, H.; Sánchez De Merás, A.; Pedersen, T. B. Reduced scaling in electronic structure calculations using Cholesky decompositions. _J. Chem. Phys._ 2003, _118_ , 9481–9484
  * Aquilante et al. 2007 Aquilante, F.; Lindh, R.; Bondo Pedersen, T. Unbiased auxiliary basis sets for accurate two-electron integral approximations. _J. Chem. Phys._ 2007, _127_ , 114107
  * Aquilante et al. 2009 Aquilante, F.; Gagliardi, L.; Pedersen, T. B.; Lindh, R. Atomic Cholesky decompositions: A route to unbiased auxiliary basis sets for density fitting approximation with tunable accuracy and efficiency. _J. Chem. Phys._ 2009, _130_ , 154107
  * Medves et al. 2022 Medves, M.; Fronzoni, G.; Stener, M. Optimization of density fitting auxiliary Slater-type basis functions for time-dependent density functional theory. _J. Comput. Chem._ 2022, _43_ , 1923–1935
  * Löwdin 1970 Löwdin, P. O. On the Nonorthogonality Problem. _Adv. Quantum Chem._ 1970, _5_ , 185–199
  * Lehtola et al. 2020 Lehtola, S.; Blockhuys, F.; Van Alsenoy, C. An overview of self-consistent field calculations within finite basis sets. _Molecules_ 2020, _25_ , 1218
  * Wirz et al. 2017 Wirz, L. N.; Reine, S. S.; Pedersen, T. B. On Resolution-of-the-Identity Electron Repulsion Integral Approximations and Variational Stability. _Journal of Chemical Theory and Computation_ 2017, _13_ , 4897–4906, PMID: 28873316
  * Langreth and Perdew 1977 Langreth, D. C.; Perdew, J. P. Exchange-correlation energy of a metallic surface: Wave-vector analysis. _Phys. Rev. B_ 1977, _15_ , 2884–2901
  * Kaltak et al. 2014 Kaltak, M.; Klimeš, J.; Kresse, G. Low scaling algorithms for the random phase approximation: Imaginary time and laplace transformations. _J. Chem. Theory Comput._ 2014, _10_ , 2498–2507
  * Kaltak et al. 2014 Kaltak, M.; Klimeš, J.; Kresse, G. Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si. _Phys. Rev. B_ 2014, _90_ , 054115
  * Helmich-Paris and Visscher 2016 Helmich-Paris, B.; Visscher, L. Improvements on the minimax algorithm for the Laplace transformation of orbital energy denominators. _J. Comput. Phys._ 2016, _321_ , 927–931
  * Förster and Visscher 2021 Förster, A.; Visscher, L. GW100: A Slater-Type Orbital Perspective. _J. Chem. Theory Comput._ 2021, _17_ , 5080–5097
  * Kozuch and Martin 2011 Kozuch, S.; Martin, J. M. DSD-PBEP86: In search of the best double-hybrid DFT with spin-component scaled MP2 and dispersion corrections. _Phys. Chem. Chem. Phys._ 2011, _13_ , 20104–20107
  * Cooper et al. 2017 Cooper, B.; Girdlestone, S.; Burovskiy, P.; Gaydadjiev, G.; Averbukh, V.; Knowles, P. J.; Luk, W. Quantum Chemistry in Dataflow: Density-Fitting MP2. _J. Chem. Theory Comput._ 2017, _13_ , 5265–5272
  * 141 Baerends, E.; Ziegler, T.; Atkins, A.; Autschbach, J.; Baseggio, O.; Bashford, D.; Bérces, A.; Bickelhaupt, F.; Bo, C.; Boerrigter, P.; Cavallo, L.; Daul, C.; Chong, D.; Chulhai, D.; Deng, L.; Dickson, R.; Dieterich, J.; Ellis, D.; van Faassen, M.; Fan, L.; Fischer, T.; Förster, A.; Guerra, C. F.; Franchini, M.; Ghysels, A.; Giammona, A.; van Gisbergen, S.; Goez, A.; Götz, A.; Groeneveld, J.; Gritsenko, O.; Grüning, M.; Gusarov, S.; Harris, F.; van den Hoek, P.; Hu, Z.; Jacob, C.; Jacobsen, H.; Jensen, L.; Joubert, L.; Kaminski, J.; van Kessel, G.; König, C.; Kootstra, F.; Kovalenko, A.; Krykunov, M.; van Lenthe, E.; McCormack, D.; Michalak, A.; Mitoraj, M.; Morton, S.; Neugebauer, J.; Nicu, V.; Noodleman, L.; Osinga, V.; Patchkovskii, S.; Pavanello, M.; Peeples, C.; Philipsen, P.; Post, D.; Pye, C.; Ramanantoanina, H.; Ramos, P.; Ravenek, W.; Rodríguez, J.; Ros, P.; Rüger, R.; Schipper, P.; Schlüns, D.; van Schoot, H.; Schreckenbach, G.; Seldenthuis, J.; Seth, M.; Snijders, J.; Solà, M.; Stener, M.; Swart, M.; Swerhone, D.; Tognetti, V.; te Velde, G.; Vernooijs, P.; Versluis, L.; Visscher, L.; Visser, O.; Wang, F.; Wesolowski, T.; van Wezenbeek, E.; Wiesenekker, G.; Wolff, S.; Woo, T.; Yakovlev, A. ADF2019.3 (locally modified), SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, https://www.scm.com
  * Dunning 1989 Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. _The Journal of Chemical Physics_ 1989, _90_ , 1007–1023
  * Pritchard et al. 2019 Pritchard, B. P.; Altarawy, D.; Didier, B.; Gibson, T. D.; Windus, T. L. New Basis Set Exchange: An Open, Up-to-Date Resource for the Molecular Sciences Community. _J. Chem. Inf. Model._ 2019, _59_ , 4814–4820
  * Van Lenthe and Baerends 2003 Van Lenthe, E.; Baerends, E. J. Optimized Slater-type basis sets for the elements 1–118. _Journal of Computational Chemistry_ 2003, _24_ , 1142–1156
  * Franchini et al. 2014 Franchini, M.; Philipsen, P. H. T.; van Lenthe, E.; Visscher, L. Accurate Coulomb Potentials for Periodic and Molecular Systems through Density Fitting. _Journal of Chemical Theory and Computation_ 2014, _10_ , 1994–2004, PMID: 26580526
  * Franchini et al. 2013 Franchini, M.; Philipsen, P. H. T.; Visscher, L. The Becke Fuzzy Cells Integration Scheme in the Amsterdam Density Functional Program Suite. _Journal of Computational Chemistry_ 2013, _34_ , 1819–1827
  * Helgaker et al. 1997 Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-set convergence of correlated calculations on water. _J. Chem. Phys._ 1997, _106_ , 9639–9646
  * Hättig 2005 Hättig, C. Optimization of auxiliary basis sets for RI-MP2 and RI-CC2 calculations: Core-valence and quintuple-$\zeta$ basis sets for H to Ar and QZVPP basis sets for Li to Kr. _Phys. Chem. Chem. Phys._ 2005, _7_ , 59–66
  * Hill and Platts 2008 Hill, J. G.; Platts, J. A. Auxiliary basis sets for density fitting-MP2 calculations: Nonrelativistic triple- all-electron correlation consistent basis sets for the 3d elements Sc-Zn. _J. Chem. Phys._ 2008, _128_ , 044104
  * 150 We have omitted the systems with index 15 and 16 which contain Iodine atoms and for which a comparison to TURBOMOLE would be difficult due to the use of pseudopotentials in ref. 50. We also do not show results for 3 of the charged systems due to convergence difficulties.
  * Petersson and Nyden 1981 Petersson, G. A.; Nyden, M. R. Interference effects in pair correlation energies: Helium L-limit energies. _J. Chem. Phys._ 1981, _75_ , 3423–3425
  * Irmler and Grüneis 2019 Irmler, A.; Grüneis, A. Particle-particle ladder based basis-set corrections applied to atoms and molecules using coupled-cluster theory. _J. Chem. Phys._ 2019, _151_ , 104107
  * Irmler et al. 2019 Irmler, A.; Gallo, A.; Hummel, F.; Grüneis, A. Duality of Ring and Ladder Diagrams and Its Importance for Many-Electron Perturbation Theories. _Phys. Rev. Lett._ 2019, _123_ , 156401
  * Petersson et al. 1988 Petersson, G. A.; Bennett, A.; Tensfeldt, T. G.; Al-Laham, M. A.; Shirley, W. A.; Mantzaris, J. A complete basis set model chemistry. I. The total energies of closed-shell atoms and hydrides of the first-row elements. _J. Chem. Phys._ 1988, _89_ , 2193–2218
  * Wu and Truhlar 2021 Wu, D.; Truhlar, D. G. How Accurate Are Approximate Density Functionals for Noncovalent Interaction of Very Large Molecular Systems? _J. Chem. Theory Comput._ 2021, _17_ , 3967–3973
  * Brémond et al. 2022 Brémond, E.; Li, H.; Sancho-García, J. C.; Adamo, C. Double Hybrids and Noncovalent Interactions: How Far Can We Go? _J. Phys. Chem. A_ 2022, _126_ , 2590–2599