Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2, and \textit{GW} with numeric atom-centered orbital basis functions

Efficient implementations of electronic structure methods are essential for first-principles modeling of molecules and solids. We here present a particularly efficient common framework for methods beyond semilocal density-functional theory, including Hartree-Fock (HF), hybrid density functionals, random-phase approximation (RPA), second-order M{\o}ller-Plesset perturbation theory (MP2), and the $GW$ method. This computational framework allows us to use compact and accurate numeric atom-centered orbitals (popular in many implementations of semilocal density-functional theory) as basis functions. The essence of our framework is to employ the"resolution of identity (RI)"technique to facilitate the treatment of both the two-electron Coulomb repulsion integrals (required in all these approaches) as well as the linear density-response function (required for RPA and $GW$). This is possible because these quantities can be expressed in terms of products of single-particle basis functions, which can in turn be expanded in a set of auxiliary basis functions (ABFs). The construction of ABFs lies at the heart of the RI technique, and here we propose a simple prescription for constructing the ABFs which can be applied regardless of whether the underlying radial functions have a specific analytical shape (e.g., Gaussian) or are numerically tabulated. We demonstrate the accuracy of our RI implementation for Gaussian and NAO basis functions, as well as the convergence behavior of our NAO basis sets for the above-mentioned methods. Benchmark results are presented for the ionization energies of 50 selected atoms and molecules from the G2 ion test set as obtained with $GW$ and MP2 self-energy methods, and the G2-I atomization energies as well as the S22 molecular interaction energies as obtained with the RPA method.


Introduction
Accurate quantum-mechanical predictions of the properties of molecules and materials (solids, surfaces, nano-structures, etc.) from first principles play an essential role in chemistry and condensed-matter research today. Of particular importance are computational approximations to the many-body Schrödinger or Dirac equations that are tractable and yet retain quantitatively reliable atomic-scale information about the system -if not for all possible materials and properties, then at least for a relevant subset. Density-functional theory (DFT) [1,2] is one such successful avenue. It maps the interacting many-body problem onto an effective single-particle one where the many-body complexity is hidden in the unknown exchange-correlation (XC) term, which has to be approximated in practice. Existing approximations of the XC term roughly fall into a hierarchial scheme [3]. Its local-density (LDA) [2] and generalized gradient approximations (GGAs) [4,5,6,7] are now well recognized workhorses with a broad application range in computational molecular and materials science. However, several qualitative failures are well known: To name but a few, certain adsorbate geometries [8], f -electron systems [9,10,11,12,13], or van der Waals interactions [14,15,16,17,18,19] are not described correctly at this level of theory. Thus, there is much ongoing work to extend the reach of density-functional theory, e.g., meta-GGAs [20,21,22], formalisms to include van der Waals interactions [23,24,25,26,27,28], hybrid functionals [29,30,31,32,33,34], or approaches based on the random-phase approximation (RPA) [35,36,37,38,39,40,41,42] that deal with the non-local correlations in a more systematical and non-empirical way.
Another avenue are the approaches of quantum chemistry, that start with Hartree-Fock theory [43,44]. These approaches offer a systematically convergable hierarchy of methods by construction, said to reach "gold standard" accuracy for many molecular systems at the level of coupled-cluster theory [45,46,47] [often, taken to include singles, doubles, and perturbative triples, CCSD(T) [48]]. CCSD(T) theory is significantly more accurate than DFT-LDA/GGA for many molecular systems but also significantly more costly (formally scaling as O(N 7 ) with system size). It has its own shortcomings as well. For instance, systematic, material-specific failures can occur in cases where the underlying Hartree-Fock solution itself is not a good reference to start with (for example, many open-shell systems), and a multireference extension of the approach [49] becomes necessary.
A third avenue for electronic structure calculations is the quantum Monte Carlo (QMC) method, in particular the diffusion QMC method [50,51]. This is a stochastic approach that deals with the many-body wavefunction directly. The diffusion QMC method can often deliver high accuracies, and provide data and insights for problems which are difficult for other approaches. Its widespread use, however, is also impeded by the rather high computational costs. Moreover, the fixed-node approximation and the underlying pseudopotential approximation are known issues which limit the practical accuracy of QMC. Regarding the computational cost, the QMC methods, the algorithm of which is intrinsically parallel, are in a better position to benefit from the development of petaflop supercomputers [52].
With the successes, but also the failures or shortcomings of the aforementioned avenues, much attention is currently devoted to the construction of further, systematic and generally applicable methods or theoretical frameworks that can offer better accuracies than conventional DFT, but have lower numerical costs and are free of the limitations of CCSD(T) and QMC. Among the various possible pathways, many-body perturbation theory (MBPT) based on an efficiently attainable and trustful electronic reference state offers such an avenue. In particular, approaches based on the RPA, which bridge the DFT and MBPT worlds [35,36,53,16], have recently enjoyed considerable attention for ground-state total-energy calculations. For electron addition and removal energies, a self-energy based approach that is consistent with the RPA total-energy treatment is Hedin's GW approximation [54]. GW is especially popular in the solid state community [55,56,57,58] and has become the method of choice for the calculation of quasiparticle band structures as measured in direct and inverse photoemission [59,60,61].
Although RPA and GW are receiving much attention in the community today, the systematic investigation of diagrammatic perturbation theory from first principles for real materials is only just beginning. Its full promise lies in the fact that it is intermediate in cost between DFT and coupled-cluster theories, and applicable in practice to molecular and condensed materials alike -including open-shell systems and metals.
Besides the more generally applicable RPA and GW approaches, another correlation method that is widely used in computational chemistry is second-order Møller-Plesset theory (MP2) [62,44], which belongs to the category of Hartree-Fock based quantum chemistry approaches mentioned above. MP2 does not reach the CCSD(T) accuracy, but its more favorable computational scaling makes it applicable to larger system sizes. In analogy to the GW self-energy, a MP2 self-energy approach [63,44] that is compatible with the MP2 total energy is possible as well. As will be demonstrated more clearly in the next section, RPA, GW , and MP2 (both total and self-energy) are related approaches both diagramatically and numerically. The development of numerical frameworks that enables the implementation of all these approaches on an equal footing with high numerical efficiency and accuracy is thus highly desirable.
In the present work, we describe the underpinnings of such a unified numerical framework that is promising to boost the efficiency for all the above-mentioned methods, by allowing for their implementations with compact and efficient NAO basis sets. Our specific implementation is based on the "Fritz Haber Institute ab initio molecular simulations" (FHI-aims) [64] program package. While we make reference to FHI-aims basis sets throughout much of this work, the numerical foundation presented here is general: applicable to any other type of atom-centered basis set. We note that many production-quality implementations of hybrid functionals, Hartree-Fock, and MBPT are based on analytical basis functions such as Gaussian-type orbitals (GTOs) or plane waves, and typically rely on pseudopotential-type approaches. In contrast, we here aim for an all-electron, full-potential treatment with NAO basis sets that does not sacrifice accuracy compared to the alternatives.
For DFT-LDA/GGA, NAOs are well established and can be found in several implementations [65,66,67,68,64,69,70]. This is however not the case for HF and the MBPT approaches we are going to address in this paper. Specifically we will present in the following i) an atom-centered resolution of identity (RI) framework analogous to what is pursued in quantum chemical methods [71,72,73,74,75,76,77,78]. This framework allows us to reduce all four-center two-electron Coulomb integrals to precomputed three-and two-center integrals. Our scheme differs from the quantum chemistry approach [77,78] in the auxiliary basis set construction, which is essential for retaining the flexibility to work with any atom-centered basis function shape, rather than being restricted to analytical shapes only.
ii) an assessment of the accuracy of the NAO basis sets used for normal LDA/GGA calculations in FHI-aims, and intended to be transferable regardless of the specific underlying materials or functionals, for Hartree-Fock, MP2, hybrid functionals, RPA, and GW .
Reference to relevant work by other groups in electronic-structure theory is made throughout this work. The present paper demonstrates our approach for molecular systems (non-periodic), and makes extensive use of established GTO basis sets for comparison and reference purposes. In addition, we provide benchmark GW vertical (geometry of the ionized molecule not relaxed) ionization energies (IEs) for a subset of the G2 ion test set [79], and benchmark binding energies (RPA) for the G2-I and the S22 molecular test set [80]. We restrict ourselves to algorithms that have standard scaling with system size [O(N 4 ) for HF, O(N 5 ) for MP2, etc]. Regarding total energy differences we restrict ourselves to a discussion of counterpoise-corrected [81] results when assessing the accuracy of MBPT methods that utilize the full (also unoccupied) spectrum.
In the following, we first recapitulate the HF, MP2, RPA, and GW methods and highlight the structural similarity and difference of the three correlation methods (section 2). We then introduce the basics of RI and the RI formulation of the above methods (section 3). Our own RI prescription and its accuracy is the subject of section 4. Section 5 demonstrates the overall accuracy of our approach with NAO basis sets for a variety of test systems for HF, MP2, hybrid density functionals, RPA, and GW . The benchmark results for a subset of the G2 and S22 molecular test sets using our approaches are presented in section 6. Finally we conclude our paper in section 7.

Many-electron Hamiltonian and many-body perturbation theory
Hartree-Fock (HF) theory, hybrid density functionals, and MBPT (MP2, RPA, and GW ) are all approximate ways to solve the interacting many-electron Hamiltonian where N e is the number of electrons that interact via the Coulomb interaction v ee ij ≡ 1/|r i − r j |, and v ext i ≡ v ext (r i ) is a local, multiplicative external potential, usually due to the nuclei. Hartree atomic units are used throughout this paper. The numerical cost for an exact solution of the Hamiltonian (1) scales exponentially with system size (number of electrons). The systems for which such a solution is possible are thus heavily restricted in size. In general, accurate approximations are needed. The most common approximations first resort to the solution of a mean-field, non-interacting Hamiltonian H 0 that yields an approximate ground-state wave function |Φ 0 : The "(0)" in E (0) 0 implies the fact that this is the ground state energy of the mean-field Hamiltonian. A suitableĤ 0 should be solvable with relative ease. |Φ 0 is a single Slater determinant formed from the lowest N e single-particle spin-orbitals determined bŷ where σ denotes the spin index andĥ 0 is the effective single-particle Hamiltonian noted in the bracket of (2). The form of (2-3) is, of course, precisely that of Kohn-Sham(KS) DFT (with a local mean-field potential v MF i ), or of HF theory and hybrid functionals (with a non-local mean-field potential v MF i ). The purpose of starting with (1-3) is to establish our notation for the following sections and to distinguish between •Ĥ (the many-electron Hamiltonian), • the mean-field HamiltonianĤ 0 , the solutions of which are many-electron wave functions given by single Slater determinants, and define an excited-state spectrum of their own (obviously not the same as that ofĤ), • the effective single-particle Hamiltonianĥ 0 , which generates the single-particle orbitals ψ nσ and orbital energies ǫ nσ .
In MBPT, one starts fromĥ 0 and the associated eigenenergies and eigenfunctions to systematically approximate the properties ofĤ, e.g., its true ground-state energy E 0 in MP2 or RPA or its singple-particle excitations in GW . The interacting many-electron HamiltonianĤ is partitioned into a mean-field HamiltonianĤ 0 as given by (2) and an interaction termĤ ′ , In the remainder of this section we collect the basic formulae that define the meanfield Hamiltonians, perturbation theory for ground state properties (MP2, RPA), and perturbation theory for excited states (electron addition and removal energies, through either GW or MP2 self-energies). From a numerical point of view, the underpinning of all these methods is the same: an efficient, accurate basis set prescription, and an efficient expansion of the two-electron Coulomb operator, which is the primary focus of this paper.

Mean-field Hamiltonians of HF or DFT
In HF theory the ground-state wave function of the Hamiltonian in (1) is approximated by a single Slater determinant |Φ 0 and E (0) 0 is obtained by a variational optimization, leading to f here denotes the HF single-particle Hamiltonian, and v h (r) is the Hartree potential, with the electron density and Σ x σ is the non-local, exact-exchange potential Equations (5)-(8) form a set of non-linear equations that have to be solved iteratively. v h (r) and Σ x σ (r, r ′ ) together yield the HF potential v HF , a special case of the meanfield potential v MF i in (2) and (4). The HF wavefunction |Φ 0 is given by the Slater determinant formed by the N e spin-orbitals ψ nσ with lowest energies ǫ nσ .
At self-consistency, the HF total energy is where the Hartree energy E h and exact-exchange energy E x are given respectively by In general, the single-particle spin-orbitals ψ nσ (r) are expanded in terms of a set of basis where c i nσ are the expansion coefficients. In terms of these basis functions, the HF effective potential can be expressed in a matrix form where in particular the exact-exchange matrix is given by In (12) D kl,σ is the density matrix and (ij|kl) is the short-hand notation of quantum chemistry for 4-center 2-electron integrals Very similar equations arise in KS-DFT with a local v MF i , or in the generalized KS scheme [82] with a fraction of Σ x σ (r, r ′ ) in the potential. In principle, the exact KS-DFT would yield the exact many-electron ground-state energy E 0 and ground-state density n 0 (r). In practice, the XC energy functional and potential have to be approximated. The effective single-particle orbitals from either HF or from approximate KS-DFT are convenient starting points for MBPT.

Perturbation theory for the many-electron ground-state energy: MP2
Assuming thatĤ ′ is benign and can be treated as a perturbation on top ofĤ 0 , the ground-state energy for the interacting system can be obtained using Rayleigh-Schrödinger perturbation theory (RSPT), or more precisely Brueckner-Goldstone perturbation theory [83,84].
According to the Goldstone theorem [84], in a diagrammatic expansion of the ground-state total energy, only the "linked" diagrams need to be taken into account. And this guarantees the the size-extensivity of the theory, i.e., the total energy scales correctly with the system size. Møller-Plesset (MP) perturbation theory is a special case of RSPT [62], where the reference HamiltonianĤ 0 the HF HamiltonianĤ HF = if i . Terminating the expansion at second order gives the MP2 theory, in which the (second-order) correlation energy is given by Here Φ k are the Slater determinants representing the excited states of H 0 = H HF , and E (0) k are the corresponding excited-state energies. H ′ is given by (4) with v MF = v HF . Among all possible excited-state configurations Φ k , only double excitations contributes in (15). This is because singly-excited Φ k do not couple to the groundstate Φ 0 (Brillouin's theorem [44] for the HF reference), whereas even higher-excited configerations (triples, quadruples, etc.) do not contribute due to the two-particle nature of the operator H ′ . As such, equation (15) can be expressed in terms of single-particle spin-orbitals, where (ma, σ|nb, σ ′ ) are 2-electron Coulomb repulsion integrals for molecular orbitals The two terms in (16) correspond to the 2nd-order Coulomb (direct) and 2nd-order exchange energy, respectively.

2.4.
Perturbation theory for the many-electron ground-state energy: RPA MP2 corresponds to the 2nd-order term in a perturbation theory where the perturbation expansion is essentially based on the bare Coulomb operator (with HF effective potential subtracted). As such, the MP2 correlation energy diverges for the homogeneous electron gas and metals with zero direct energy gap. To overcome this problem in the framework of MBPT it is essential to sum up the diverging terms in the perturbation series to infinite order. One such example, which has gained considerable popularity recently [38,85,86,87,88,89,90,91,92,39,93,94,95,96,97,98], is the RPA [35,36,99,53,16], that in the context of MBPT corresponds to an infinite summation of "ring" diagrams. The choice of the non-interacting reference HamiltonianĤ 0 can be, e.g., HF or DFT with any desired XC functional. Apart from the diagrammatic representation, RPA can also be formulated in other ways, e.g., as the simplest time-dependent Hartree approximation in the context of the adiabatic-connection fluctuation-dissipation theorem [99,53,16], or as a subclass of terms in coupled-cluster theory with double excitations [87]. In the context of DFT [53], RPA calculations can be performed self-consistently by means of the optimized effective potential approach [100,101,102,103].
In close-packed notation, the RPA correlation energy (cRPA) can be expressed in terms of the RPA dielectric function ε, or alternatively the non-interacting density response function χ 0 on the imaginary frequency axis The real-space (Adler-Wiser [104,105]) representation of χ 0 reads where c.c. denotes "complex conjugate", and ψ n (r) and ǫ n are single-particle orbitals and orbital energies as implied by (3). The RPA dielectric function ε is linked to χ 0 through Using (17) and (19), the lowest-order term in (18) can be expressed as Higher-order terms in (18) follow analogously. There is thus a straightforward path to compute χ 0 , and hence the RPA correlation energy, once the selected non-interacting reference state is solved. In practice the RPA correlation energy is always combined with the exact-exchange energy (EX) in (9), henceforth denoted (EX+cRPA), but evaluated with the same single-particle orbitals as used in E RPA c . The choice of input orbitals will in the following be marked by (EX+cRPA)@MF, where MF specifies the mean-field approach used to compute the single-particle orbitals. The application of EX+cRPA to various realistic systems, as well as the development of schemes that go beyond simple EX+cRPA, is an active field [38,85,86,106,107,87,88,91,108,92,39,89,93,94,95,96,97,98,109,40]. For instance, when the so-called single-excitation and secondorder screened exchange contributions [110,90] are added to EX+cRPA, the resulting accuracy is impressive [40,41].

Perturbation theory for electron addition or removal energies: GW or MP2
Ground-state energies aside, one is often interested in the properties of electronically excited states. Part of this information is in principle accessible by taking the difference of the total energies of N-electron system and N ± 1-electron systems using approaches discussed above. In practice, this approach mainly works for computing core-level excitations and/or the first ionization energy and electron affinity for (small) finite systems. In contrast, Green function techniques are more convenient and powerful for dealing with electronic excitations in general. The basic theory of Green functions is well documented in textbooks [111,112]. Here we collect the contextual equations, based on which practical approximations can be introduced.
The single-particle Green function of an interacting many-electron system is defined as G(r, t; r ′ , t ′ ) = −i N|Tψ(r, t)ψ † (r ′ , t ′ )|N (22) where |N = |Ψ 0 (r 1 , . . . , r N ) denotes the interacting ground-state wave function of an N-electron system (solution to (1)).ψ(r, t) andψ † (r ′ , t ′ ) are field operators in the Heisenberg picture that annihilate and create an electron at space-time point (r, t) and (r ′ , t ′ ), respectively.T is the time-ordering operator. The Green function G(r, t; r ′ , t ′ ) measures the probability amplitude of a hole created at (r, t) propagating to (r ′ , t ′ ) for t < t ′ , or an electron added at (r ′ , t ′ ) propagating to (r, t) for t > t ′ . The poles of its Fourier transform, G(r, r ′ , ω), correspond to the single-particle excitation energies as measured for example in direct and inverse photoemission experiments. For Hamiltonians with a time-independent external potential, as considered in this paper, the Green function depends only on the difference between t and t ′ , G(r, t; r ′ , t ′ ) = G(r, r ′ ; t−t ′ ). Its Fourier transform gives the frequency-dependent Green function G(r, r ′ ; ω) that satisfies the Dyson equation, Here v h (r) is the electrostatic Hartree potential defined in (6) and Σ(r, r ′′ , ω) is the dynamical, non-local, complex self-energy that contains all the many-body XC effects. G(r, r ′ ; ω) and Σ(r, r ′ ; ω) for many-body interacting systems can in principle be obtained using diagramatic Feynman-Dyson perturbation theory. The perturbation series is built on a non-interacting Green function G 0 (r, r ′ , ω) that corresponds to the non-interacting single-particle Hamiltonianĥ 0 . With single-particle orbitals ψ nσ (r) and orbital energies ǫ nσ ofĥ 0 , one has where ǫ F is the Fermi energy, and η a positive infinitesimal. For a given G 0 , corrections to the single-particle excitation energies can be computed from approximate perturbative expansions of Σ(r, r ′ ; ω). Examples are the GW method and the 2nd-order approximations discussed below.
In the GW approximation proposed by Hedin [54], the self-energy assumes the form Here W is the screened Coulomb potential at the RPA level with the dynamical dielectric function ε as defined in (19) and (20), but on the real frequency axis with iω replaced by ω + iη (η → 0 + ). In practice, one-shot perturbative GW calculations (often referred to as G 0 W 0 ) based on a fixed, non-interacting reference state (DFT with popular functionals or HF) are often performed. With χ 0 computed from the non-interacting Green function in (24), the ensuing W can be expanded in powers of The GW approximation can thus be regarded as an infinite series in the bare Coulomb potential v, or alternatively, as is obvious from (25), as first-order perturbation in terms of the screened Coulomb potential W . Once the G 0 W 0 self-energy is obtained from (25), the corrections to the singleparticle orbital energies are given by where the XC part v xc of the reference mean-field potential has to be subtracted. 28 approximates the quasiparticle wavefunctions with the single-particle orbitals of the reference state. This is often justified, but may break down in certain cases [113,114,115,116,117,118].
Rewriting the diagonal elements of the G 0 W 0 self-energy on the imaginary frequency axis in terms of four-center Coulomb integrals, we obtain where G 0 mσ (iω) = 1/(iω + ǫ F − ǫ mσ ), and This expression can be analytically continued to the real-frequency axis using Padé approximation or a two-pole model [119]. The G 0 W 0 approach, as a highly popular choice for quasiparticle excitation calculations, in particular for solids [59,120,60,61,121], is akin to the RPA approach for computing the ground-state correlation energy. Similarly, one can also introduce a 2ndorder perturbation theory for the self-energy [44] consistent with the MP2 correlation energy. Based on the HF noninteracting Green function [equation (24) with HF orbitals], the 2nd-order self-energy can be expressed as Here for notational simplicity we have used 1 = (r 1 , t 1 ), and v(1, 2) = v(r 1 − r 2 )δ(t 1 , t 2 ). After integrating out the internal (spatial and time) coordinates in (31), the final expression for Σ (2) (in frequency domain) within the HF molecular orbital basis is given by where θ(x) is the Heaviside step function and η → 0 + . Again, the two lines in (32) correspond to correlations arising from the 2nd-order direct and 2nd-order exchange interaction, respectively. With this self-energy single-particle excitation energies -here denoted MP2 quasiparticle energies -can be obtained by adding a correction to the HF orbital energies Similar to the MP2 correlation energy, the quality of the MP2 self-energy as given by (32) will deteriorate as the single-particle energy gap shrinks. By means of the Dyson equation G = G 0 + G 0 ΣG self-consistent Green functions and self-energies could be obtained. The advantage is that the self-consistent Green function will then be independent of G 0 and satisfies particle number, momentum and energy conservation [122,123]. For an in-depth discussion, we refer to [124,125] and references therein.
3. Resolution of identity for HF, MP2, RPA, and GW

Background
In this section we present the basic resolution of identity (RI) formalism: the auxiliary basis sets, different variants of RI, and the working equations for HF, hybrid density functionals, MP2, RPA, and GW . Similar accounts have been given in the literature on one or another of the above methods, and we encourage readers to consult them [76,77,78,126,127,128,129]. Our aim here is to lay out a complete description of all necessary specifics in one consistent notation that will allow us to present our own developments (next sections) in a self-contained way.
The common ingredient of all techniques introduced in section 2 (Hartree-Fock, hybrid functionals, MP2, RPA, GW ) are the four-orbital Coulomb integrals where (ij|kl) are the two-electron integrals in a basis set representation as defined in (14), and c i mσ are the eigen-coefficients for the molecular orbitals. Computing (and possibly storing) these 4-center 2-electron integrals can be a major bottleneck for approaches beyond LDA and GGAs.
For analytical GTOs algorithms have been developed to handle (mn, σ|ab, σ ′ ) integrals efficiently and on-the-fly [130,131]. More general NAOs, however, are not amenable to such algorithms. In the context of HF, we note that RI is not the only technique available to deal with the four-center integrals: Making use of the translational properties of spherical harmonics, Talman and others [132,133,134,135] have developed techniques based on multipole expansions of basis functions. Multi-center NAO integrals can then be treated partially analytically. Alternatively, efficient Poisson solvers [136] have recently been used to enable direct NAO HF calculations through four-center integrals for simple systems [133,136]. Finally, a yet different numerical route (based on expanding orbital products directly) has been adopted along the line of time-dependent DFT in the linear-response framework [137]. RI is, however, most successful at reducing the computational load compared to direct four-center integral based methods, most prominently for MP2 in quantum chemistry [77], which is why we pursue this route here.

Auxiliary basis
The resolution of identity (RI) or synonymously density fitting technique [71,72,73,74,75,76,77,78] amounts to representing pair products of atomic basis functions ϕ i (r)ϕ j (r) in terms of auxiliary basis functions (ABFs), µ = 1, 2, . . . , N aux labels the auxiliary basis functions {P µ }, C µ ij are the expansion coefficents, and ρ ij (r) andρ ij (r) here denote pair products of basis functions and their approximate expansion in ABFs. The evaluation of the 4-center integrals in (14) then reduces to To determine the expansion coefficients C µ ij , three-center integrals involving the ABFs and the pair products of the NAOs are required. Thus the expensive (both in time and memory, if there is a need to pre-compute numerical matrix elements) 4-center integrals reduce to the much cheaper 3-center and 2-center ones in RI. The key reason for the success of RI lies in the fact that the set of all possible pair products {ϕ i (r) · ϕ j (r)}, as a set of basis functions in three-dimensional function space, is heavily linearly dependent. Their number scales quadratically with system size, while a non-redundant basis set that expands the same three-dimensional space should scale linearly with system size. For example, the non-interacting response function χ 0 in (19), as well as the screened Coulomb interaction W in (26), is written in terms of orbital pair products, and hence can be represented in terms of the ABFs. As will be shown below, this naturally leads to a RI implementation for RPA and GW .
Next we will present RI formulations for all pertinent methods in this paper before presenting our specific choice for the ABFs subsequently.

Metric and variational principle in RI
For a given set of ABFs {P µ (r)}, the way to determine the expansion coefficients C µ ij is not unique. Different variational procedures give rise to different versions of RI and different working equations for computing the C µ ij [71,72,73,74,75,76,77,78]. The expansion error of basis products in terms of the ABFs [equation (35)] is One choice for the construction of the expansion coefficients C µ ij is to minimize the residual δρ ij (r). A simple least-square fit amounts to minimizing the norm of the residual, |δρ ij (r)| 2 dr and yields where and Combining (36) and (39), one arrives at the following approximation to the four-center Coulomb integrals, In the literature equation (42) is therefore referred to as the "SVS" version [75] of RI ("RI-SVS" in the following) because of the appearance of the inverse S before and after the V matrix in (42). A better criterion for obtaining C µ ij is to directly minimize the RI error of the 4-center integrals themselves, As shown by Whitten [72] δI ij,kl has an upper bound: The minimization of δI ij,kl can thus be achieved by independently minimizing δU ij = (δρ ij |δρ ij ) and δU kl = (δρ kl |δρ kl )-the self-repulsion of the basis pair density residuals. Minimizing δU ij with respect to C µ ij leads to [72,73,74,75] where and the Coulomb matrix V is defined in (37). Combining (36) and (45) one obtains the following decomposition of the 4-center ERIs Equation (47) is based on the global Coulomb metric and corresponds to the "RI-V" method [75] mentioned earlier. "RI-V" has long been known to be superior to "RI-SVS" in the quantum chemistry community [72,73,74,75]. This can be most easily seen by inspecting the error introduced in the self-Coulomb-repulsion of the NAO pairs, In the RI-V approximation, the first term in the above equation vanishes, and the nonzero contribution comes only from the second order of δρ ij . This can be readily verified as follows, where (35) and particularly (45) have been used. In contrast, in RI-SVS the term linear in δρ is non-zero and represents the dominating contribution to the total error. Our preferred flavor of RI is therefore RI-V, based on the Coulomb metric [72,73,74,75,76,77,78], on which all working equations for HF and other approaches presented further down in this section are based. Before proceeding we reiterate that RI-V continues to be the de facto standard in quantum chemical calculations, due to its well-established accuracy and reliability [77,78].
That said, the long-range nature of the Coulomb interaction does present a bottleneck for implementations that scale better than the textbook standard (e.g., better than O(N 4 ) for Hartree-Fock). In order to avoid delocalizing each localized two-center basis function product entirely across the system through C µ ij , more localized approaches would be desirable. Research into better-scaling RI expansions that retain at least most of the accuracy of the Coulomb metric is thus an active field, for example by Cholesky decomposition techniques or an explicitly local treatment of the expansion of each product (see, e.g., [138,139,140,141,142,143,144,145] for details). In fact, a promising Coulomb-metric based, yet localized, variant of RI has been implemented in FHI-aims. In this approach products of orbital basis functions are only expanded into auxiliary basis functions centered at the two atoms at which the orbital basis functions are centered, but the appropriate RI sub-matrices are still treated by the Coulomb metric [146]. As expected, the error cancellation in this approach is not as good as that in full RI-V, but-for Hartree-Fock and hybrid functionals-certainly more than an order of magnitude better than in RI-SVS, creating a competitive alternative for cases where RI-V is prohibitive. More details would go beyond the scope of this paper and will be presented in a forthcoming publication [147].

HF and hybrid functionals
The key quantity for HF and hybrid functionals is the exact-exchange matrix -the representation of the non-local exact-exchange potential [equation (8)] in terms of basis functions as given in 12). Its RI-V expansion follows by inserting (47) into (12): Σ x ij,σ must thus be recomputed for each iteration within the self-consistent field (scf) loop. The required floating point operations scale as N 3 b · N aux . The transformation matrix M µ ik , on the other hand, is constructed only once (prior to the scf loop), requiring N 2 b · N 2 aux operations. The numerical efficiency can be further improved by inserting the expression for the density matrix (13) into (49): The formal scaling in (51) is now N occ · N b 2 · N aux , with N occ being the number of occupied orbitals, and thus improved by a factor N b /N occ (typically 5 to 10). Once the exact-exchange matrix is obtained, the exact-exchange energy follows through For a variety of physical problems, combinding HF exchange with semi-local exchange and correlation of the GGA type gives much better results than with pure HF or pure GGAs [29]. Various flavors of these so-called hybrid functionals exist in the literature. The simplest one-parameter functionals are of the following form In the PBE0 hybrid functional [30], the GGA is taken to be PBE, and the mixing parameter α is set to 1/4. Naturally, the computational cost of hybrid functionals is dominated by the HF exchange. Once HF exchange is implemented, it is straightforward to also perform hybrid functional calculations.

MP2 (total-energy correction and self-energy)
To compute the MP2 correlation energy in (16) and the MP2 self-energy in (32) using the RI technique, the MO-based 4-orbital 2-electron Coulomb integrals are decomposed as follows The 3-orbital integrals can be evaluated by following (34), (47), and (50). Plugging (54) into (16), one obtains the RI-V version of the MP2 correlation energy In practice, one first transforms the atomic orbital-based integrals M µ ij to MO-based ones O µ ma,σ . The transformation scales formally as N occ · N 2 b · N aux and the summation in (56) Like in Hartree-Fock, the scaling exponent O(N 5 ) is therefore not reduced by RI-V. However, the prefactor in RI-MP2 is one to two orders of magnitude smaller than in full MP2. The computation of the MP2 self-energy at each frequency point proceeds analogously to that of the correlation energy. In our implementaton, we first calculate the MP2 self-energy on the imaginary frequency axis, Σ (2) mn,iσ , and then continue analytically to the real frequency axis using either a "two-pole" model [119], or the Padé approximation. Both approaches have been implemented in FHI-aims and can used to cross-check each other to guarantee the reliability of the final results.

RPA and GW
To derive the working equations for the RPA correlation energy and the GW self-energy in the RI approximation, it is illuminating to first consider the RI-decomposition of Tr[χ 0 (iω)v]. Combining (19) and (54) we obtain where O µ ma,σ is given by (55). Next we introduce an auxiliary quantity Π(iω): which allows us to write Thus the matrix Π can be regarded as the matrix representation of the composite quantity v 1/2 χ 0 (iω)v 1/2 using the ABFs. It is then easy to see that the RPA correlation energy (18) can be computed as using the general property Tr[ln(A)] = ln [det(A)] for any matrix A. This is very convenient since all matrix operations in 60 occur within the compressed space of ABFs and the computational effort is therefore significantly reduced.
In practice, we first construct the auxiliary quantity Π(iω) [equation (58)] on a suitable imaginary frequency grid where we use a modified Gauss-Legendre grid (see appendix 7 for further details) with typically 20-40 frequency points. For fixed frequency grid size, the number of required operations is proportional to N occ · N unocc · N 2 aux (N unocc is the number of unoccupied orbitals using the full spectrum of our Hamiltonian matrix). The next step is to compute the determinant of the matrix 1 − Π(iω) as well as the trace of Π(iω). What remains is a simple integration over the imaginary frequency axis. Thus our RI-RPA implementation is dominated by the step in (58) that has a formal O(N 4 ) scaling. An O(N 4 )-scaling algorithm of RI-RPA was recently derived from a different perspective [129], based on the plasmonic formulation of RPA correlation energy [107] and a transformation analogous to the Casimir-Polder integral [148]. An even better scaling can be achieved by taking advantage of the sparsity of the matrices involved [145].
Finally we come to the RI-V formalism for GW . To make (30) tractable, we expand the screened Coulomb interaction W (iω) in terms of the ABFs. Using (27) and To apply the RI-decomposition to (30), we expand where (10) and (35) are used. Combing (30), (50), (55), and 61) gives Inserting (62) into (29), one finally arrives at the RI-version of the G 0 W 0 self-energy As stated above for the MP2 self-energy, the expression is analytically continued to the real-frequency axis before the quasiparticle energies are computed by means of (28).

Orbital basis set definitions
For the practical implementation, all the aforementioned objects (wave functions, effective single-particle orbitals, Green function, response function, screened Coulomb interaction etc.) are expanded either in a single-particle basis set or an auxiliary basis set. We first summarize the nomenclature used for our NAO basis sets [64] before defining a suitable auxiliary basis prescription for RI in the next subsections. NAO basis sets {ϕ i (r)} to expand the single-particle spin orbitals ψ nσ (r) [equation (10)] are of the general form u s(a)lκ is a radial function centered at atom a, and Y lm (r a ) is a spherical harmonic. The index s(a) denotes the element species s for an atom a, and κ enumerates the different radial functions for a given species s and an angular momentum l. The unit vector r a = (r − R a )/|r − R a | refers to the position R a of atom a. The basis index i thus combines a, κ, l, and m. For numerical convenience, and without losing generality, we use real-valued basis functions, meaning that the Y lm (Ω) denote the real (for m = 0, · · · , l) and the imaginary part (for m = −l, · · · , −1) of complex spherical harmonics. For NAOs, u s(a)lκ (r) need not adhere to any particular analytic shape, but are tabulated functions (in practice, tabulated on a dense logarithmic grid and evaluated in between by cubic splines). Of course, Gaussian, Slater-type, or even muffin-tin orbital basis sets are all special cases of the generic shape (64). All algorithms in this paper could be used for them. In fact, we employ the Dunning GTO basis sets (see [149,150] and references therein) for comparison throughout this work.
Our own implementation, FHI-aims, [64] provides hierarchical sets of all-electron NAO basis functions. The hierarchy starts from the minimal basis composed of the radial functions for all core and valence electrons of the free atoms. Additional groups of basis functions, which we call tiers (quality levels) can be added for increasing accuracy (for brevity, the notation is minimal, tier 1, tier 2, etc.). Each higher level includes the lower level. In practice, this hierarchy defines a recipe for systematic, variational convergence down to meV/atom accuracy for total energies in LDA and GGA calculations. The minimal basis (atomic core and valence radial functions) is different for different functionals, but one could as well use, e.g., LDA minimal basis functions for calculations using other functionals in cases their minimal basis functions are not available. We discuss this possibility for HF below.
To give one specific example, consider the nitrogen atom (this case and more are spelled out in detail in Table 1 of Ref. [64]). There are 5 minimal basis functions, of 1s, 2s, and 2p orbital character, respectively. In a shorthand notation, we denote the number of radial functions for given angular momenta as (2s1p) (two s-type radial functions, one p-type radial function). At the tier 1 basis level, one s, p, and d function is added to give a total of 14 basis functions (3s2p1d). There are 39 basis functions (4s3p2d1f 1g) in tier 2, 55 (5s4p3d2f 1g) in tier 3, and 80 (6s5p4d3f 2g) in tier 4.

Construction of the auxiliary basis
In the past, different communities have adopted different strategies for building auxiliary basis sets. For the GTO-based RI-MP2 method [77], which is widely used in the quantum chemistry community, a variational procedure has been used to generate optimal gaussian-type atom-centered ABF sets. In the condensed matter community a so-called "product basis" has been employed in the context of all-electron GW implementations based on the linearized muffin-tin orbital (LMTO) and/or augmented plane-wave (LAPW) method [151] to represent the response function and the Coulomb potential within the muffin-tin spheres [126,152,127]. Finally, it is even possible to generate ABFs only implicitly, by identifying the "dominant directions" in the orbital product space through singular value decomposition (SVD) [144,145,153]. As will be illustrated below, our procedure to construct the ABFs combines features from both communities. Formally it is similarly to the "product basis" construction in the GW community, but instead of the simple overlap metric the Coulomb metric is used to remove the linear dependence of the "products" of the single-particle basis functions, and to build all the matrices required for the electronic structure schemes in this paper.
Our procedure employs numeric atom-centered ABFs whereby the infrastructure that is already available to treat the NAO orbital basis sets can be utilized in many respects. Specifically the ABFs are chosen as just like for the one-particle NAO basis functions in (64), but of course with different radial functions. To distinguish the auxiliary basis functions from the NAO basis functions we denote the radial functions of the ABFs as ξ s(a)lκ . The auxiliary basis should primarily expand products of basis functions centered on the same atom exactly, but at the same time be sufficiently flexible to expand all other two-center basis function products with a negligible error. In contrast to the ABFs used in the GTO-based RI-MP2 method [77], in our case the construction of auxiliary basis functions follows from the definition of the orbital basis set. At each level of NAO basis, one can generate a corresponding ABF set, denoted as aux min, aux tier 1, aux tier 2, etc. We achieve this objective as follows: (i) For each atomic species (element) s, and for each l below a limit l max s , we form all possible "on-site" pair products of atomic radial functions {ξ slκ (r) = u sk 1 l 1 (r)u sk 2 l 2 (r)}. The allowed values of l are given by the possible multiples of the spherical harmonics associated with the orbital basis functions corresponding to u sk 1 l 1 and u sk 2 l 2 , i.e., |l 1 − l 2 | ≤ l ≤ |l 1 + l 2 |.
(ii) Even for relatively small orbital basis sets, the number of resulting auxiliary radial functions {ξ slκ (r)} is large. They are non-orthogonal and heavily linear dependent. We can thus use a Gram-Schmidt like procedure (separately for each s and l) to keep only radial function components ξ slκ (r) that are not essentially represented by others, with a threshold for the remaining norm ε orth , below which a given radial function can be filtered out. In doing so the Coulomb metric is used in the orthogonalization procedure. The result is a much smaller set of linearly independent, orthonormalized radial functions {ξ slκ (r)} that expand the required function space.
(iii) The radial functions {ξ s(a),lκ } are multiplied with the spherical harmonics Y lm (r a ) as in (65).
(iv) The resulting {P µ (r)} are orthonormal if they are centered on the same atom, but not if situated on different atoms. Since we use large ABF sets, linear dependencies could also arise between different atomic centers, allowing us to further reduce the ABF space through SVD of the applied metric (S in the case of RI-SVS, V in the case of RI-V). For the molecule-wide SVD we use a second threshold ε svd , which is not the same as the on-site Gram-Schmidt threshold ε orth .
For a given set of NAOs, the number of the corresponding ABFs depends on the angular momentum limit l max s in step 1 and the Gram-Schmidt orthonormalization threshold ε orth , and to a small extent on ε svd . For RI-V, as documented in the literature [78] and demonstrated later in this work (section 4.5), it is sufficient to keep l max s just one higher than the highest angular momentum of the one-electron NAOs. Usually ε orth = 10 −2 or 10 −3 suffices for calculations of energy differences. Nevertheless both l max s and ε orth can be treated as explicit convergence parameters if needed. In practice, we keep ε svd as small as possible, typically 10 −4 or 10 −5 , only large enough to guarantee the absence of numerical instabilities through an ill-conditioned auxiliary basis. The resulting auxiliary basis size is typically 3-6 times that of the NAO basis. This is still a considerable size and could be reduced by introducing optimized ABF sets as is sometimes done for GTOs. On the other hand, it is the size and quality of our auxiliary basis that guarantees low expansion errors for RI-V, as we will show in our benchmark calculations below. We therefore prefer to keep the safety margins of our ABFs to minimize the expansion errors, bearing in mind that the regular orbital basis introduces expansion errors that are always present.

Numerical integral evaluation
With a prescription to construct ABFs at hand, we need to compute the overlap integrals C µ ij , defined in (39) for RI-SVS or in (45) for RI-V, respectively. We also need their Coulomb matrix given by (37) in general, and additionally the "normal" overlap matrix S µν given by (41) in RI-SVS. Having efficient algorithms for these tasks is enormously important, but since many pieces of our eventual implementation exist in the literature, we here only give a brief summary and refer to separate appendices for details.
Since our auxiliary basis set {P µ } is atom-centered, we obtain the Coulomb potential Q µ (r) of each P µ (r) by a one-dimensional integration for a single multipole (7). The required three-center integrals are carried out by standard overlapping atom-centered grids as used in many quantumchemical codes for the exchange-correlation matrix in DFT [154,65,64,155], see 7. The same strategy works for two-center integrals (µ|ν) = P µ (r)Q ν (r)dr.
As an alternative, we have also implemented two-center integrals following the ideas developed by Talman [133,134], which are described in 7 and 7. In summary, we thus have accurate matrix elements at hand that are used for the remainder of this work.

Accuracy of the auxiliary basis: expansion of a single product
In this section we examine the quality of our prescription for generating the ABFs as described in section 4.2. Our procedure guarantees that the ABFs accurately represent the "on-site" products of the NAO basis pairs by construction, but it is not a priori clear how the "off-site" pairs are represented. The purpose of this section is to demonstrate the quality of our ABFs to represent the "off-site" pairs. In the left panels of figure 1 we plot ρ 2s-2px (r) for a simple N 2 molecule at the equilibrium bonding distance (d = 1.1Å) -the product of the atomic 2s function from the left atom and the atomic 2p x function from the right atom. We compare this directly taken product to its ABF expansions, both in RI-SVS [equation (39)] and in RI-V [equation (45)]. The particular product ρ 2s-2px (r) is part of the minimal basis of free-atom like valence radial functions. As we increase the orbital basis set (adding tier 1, tier 2, etc.), the exact product remains the same, whereas its ABF expansion will successively improve, since the auxiliary basis set is implicitly defined through the underlying orbital basis. Three different levels of ABF sets are shown (aux min, aux tier 1, aux tier 2 from the top to bottom panels). The onsite threshold ε orth is set to 10 −2 , and the global SVD threshold ε svd is set to 10 −4 , yielding 28, 133, and 355 ABFs, respectively. In the right panels of figure 1, the corresponding δρ 2s-2px (r)the deviations of the ABF expansions from the reference curve -are plotted. One can clearly see two trends: First, the quality of the ABF expansion improves as the number of ABFs increases. Second, at the same level of ABF, the absolute deviation of the RI-V expansion is larger than the RI-SVS expansion. This is an expected behaviour for the simple pair product: RI-V is designed to minimize the error of the Coulomb integral, see section 3.3. In either method, the remaining expansion errors are centered around the nucleus, leading to a relatively small error in overall energies (in the 3-dimensional integrations, the integration weight r 2 dr is small).  Table 1 gives the errors of the Coulomb repulsion δI ij,ij of the 2s-2p x NAO basis pair under the RI approximation for the three ABF basis sets of figure 1. Table 1 also includes the influence of the threshold parameters ε orth and ε svd , separate for RI-SVS (top half) and RI-V (bottom half). The error diminishes quickly with increasing ABF basis size, and is 2-3 orders of magnitude smaller in RI-V than in RI-SVS. By decreasing ε orth , the number of ABFs at each level increases, improving particularly the accuracies at the aux min and aux tier 1 levels. The global SVD threshold ε svd comes into play only for the larger basis sets. In general, both control parameters have a much bigger effect on RI-SVS than RI-V, underscoring the desired variational properties of RI-V [71,72,73,74,75].

Accuracy of the auxiliary basis: energies and thresholds
Next we turn to the accuracy of our auxiliary basis prescription for actual HF and MP2 (total and binding energy) calculations. For this purpose, we employ all-electron GTO basis sets, when possible, to be able to refer to completely independent and accurate implementations from quantum chemistry without invoking the RI approximation (referred to as "RI-free" in the following). Our specific choice here is the GTO-based NWChem [156] code package, where "RI-free" results can be obtained using traditional methods of quantum chemistry. In the following we compare our RI-based HF and MP2 results with their "RI-free" counterparts produced by NWChem, in order to benchmark the accuracy of the RI implementation in FHI-aims. All results presented in this section Table 1. The errors δI 2s-2px,2s-2px (in eV) introduced in the RI-SVS and RI-V for calculating the self-repulsion of NAO pair products (ρ 2s-2px |ρ 2s-2px ) for N 2 (d = 1.1Å) at three levels of ABF basis sets. The number of ABFs that survive the SVD is also shown. correspond to the cc-pVQZ basis set of Dunning et al. [157,149,150] The convergence behaviour with respect to NAO / GTO single-particle basis set size will be the topic of next section. We first check the quality of the ABF prescription for light (first and secondrow) elements. We again choose N 2 as a first illustrating example. Table 2 presents RI-HF total and binding energy errors for the N 2 molecule at bond length d=1.1Å. The reference numbers given at the bottom of the table are from "RI-free" NWChem calculations. All other numbers were obtained with FHI-aims and the ABF prescription of section 4.2. Total and binding energy errors are given for several different choices Table 2. The deviation of HF total (∆E tot ) and binding energies (∆E b ) (in meV) from NWChem reference values for N 2 at bond length d =1.1Å. RI-V and RI-SVS calculations are done using the Dunning cc-pVQZ basis set. [157,149] The reference E tot and E b values (in eV) from NWChem calculations are shown at the bottom. The binding energies are BSSE-corrected using the counterpoise method. For N 2 at equilibrium bonding distance, the salient results can be summarized as follows. First, we see that RI-V with our ABF prescription implies total energy errors for Hartree-Fock of only ∼0.1-0.2 meV, while the corresponding RI-SVS errors are much larger. Both methods can, however, be adjusted to yield sub-meV binding energy errors, with those from RI-V being essentially zero. This is consistent with our observations for the error in the self-repulsion integrals in section 4.4. Since RI-V performs much better than RI-SVS, we will only report RI-V results for the remainder of this paper.
The excellent quality of our ABFs is not restricted to the equilibrium region of N 2 . In the top panel of figure 2 the restricted HF total energies are plotted for a range of bonding distances. The RI-V numbers are in very good agreement with the reference throughout. For greater clarity, the total-energy deviation of the RI-V result from the reference is plotted in the middle panel of figure 2. One can see that the totalenergy errors are in general quite small, but the actual sign and magnitude of the errors vary as a function of bond length. The deviation is shown for two different choices of the ABFs, (ε orth =10 −2 , ε svd =10 −4 ) (standard accuracy) and (ε orth = 10 −3 , ε svd =10 −5 ) (somewhat tighter accuracy). It is evident that the tighter settings produce a smoother total energy error at the sub-meV level, but, strikingly, there is no meaningful difference for the counterpoise corrected binding energy of N 2 (bottom panel of figure 2).
Next we demonstrate the quality of our ABFs for a set of molecules consisting of first and second-row elements. In figure 3 and 4 the non-relativistic RI-V HF and    Figure 3. (color online) Deviations of RI-V HF total energies (red circles) and atomization energies (blue squares) from the corresponding reference values for 20 small molecules. The three panels illustrate the dependence of the RI errors on the truncation parameters ε orth , ε svd , and the highest ABF angular momentum l ABF-max (l AO-max denotes the highest angular momentum of the single-particle atomic orbitals). Experimental equilibrium geometries and the gaussian cc-pVQZ basis are used.    MP2 total energy errors and atomization energy errors for 20 molecules are shown. In all cases the total energy error is below 1 meV /atom, demonstrating that the meVaccuracy in total energy can routinely be achieved for the RI-V approximation with our ABFs for light elements. In addition, it is evident that varying the auxiliary basis convergence settings has essentially no influence on the low overall residual error, which is attributed to other small numerical differences between two completely different codes (analytical integrations in NWChem vs. numerical integrations in FHI-aims, for example). Furthermore, it is also clear that it is sufficient to choose the highest angular momentum in the ABF construction (l ABF-max ) to be just 1 higher than that of the single-particle atomic orbitals (l AO-max ).
Having established the quality of our ABFs for light elements, we now proceed to check their performance for the heavier elements where some noteworthy feature is emerging. In figure 5 we plot the errors in the RI-V HF total energies and binding energies for Cu 2 as a function of the bond length. Again the cc-pVQZ basis and the NWChem reference values are used here. Using the thresholding parameters (ε orth =10 −2 , ε svd =10 −4 ), the RI-V HF total energy error can be as large as ∼ 15 meV /atom for copper, in contrast to the < 1 meV /atom total energy accuracy for light elements. This is because for Cu, deep core electrons are present and the absolute total-energy scale is 1-2 orders of magnitude larger than that of light elements. The residual basis components eliminated in the on-site Gram-Schmidt orthonormalization procedure (see section 4.2) thus give bigger contributions to the total energy on an absolute scale (although not on a relative scale). Indeed by decreasing ε orth the total-energy error gets increasingly smaller, and 1-1.5 meV /atom total-energy accuracy can be achieved at ε orth = 10 −4 , , l ABF-max =l AO-max +1 as demonstrated in the upper panel of figure 5. In contrast, similar to the case of light elements, the errors in the BSSE-corrected binding energies are significantly below 0.1 meV /atom along a large range of bonding distances, regardless of the choice of thresholding parameters. And also increasing the highest angular momentum for ABFs beyond l AO-max + 1 does not give noticeable improvements. Finally we look at the quality of our ABFs for even heavier elements -the Au dimer. Due to the strongly localized core states in Au, all-electron GTO basis sets that are converged to the same level of accuracy as for N and Cu above are, to our knowledge, not available for Au. Furthermore, relativity is no longer negligible and must at least be treated at a scalar relativistic level. However, different flavors of relativistic implementations can differ heavily in their absolute total energy (even if all chemically relevant energy differences are the same). An independent reference for all-electron GTO-based HF total energy with the same relativistic treatments available in FHI-aims is not readily obtainable. Under such circumstances, we demonstrate here the total energy convergence with respect to our own set of thresholding parameters.
In figure 6 we plot the deviations of the RI-V HF total and binding energies for Au 2 obtained with FHI-aims using NAO tier 2 basis with somewhat less tight thresholding parameters from those obtained with a very tight threshold setting (ε orth = 10 −5 , ε svd = 10 −5 ). The relativistic effect is treated using the scaled zeroth-order regular approximation (ZORA) [158] (see section 4.5), but this detail is not really important for the discussion here. From figure 6 one can see that with thresholding parameters that are sufficient for light elements (ε orth =10 −2 , ε svd =10 −4 ), the error in the total energy is even bigger -one order of magnitude larger than in the case of Cu 2 . However, by  going to tighter and tighter on-site ABF thresholding parameter ε orth , the total-energy error can again to be reduced to the meV /atom level. Similar to the Cu 2 case, the accuracy in the BSSE-corrected binding energy is still extraordinarily good, well below 0.1 meV regardless of the choice of thresholding parameters. The counterpoise correction can thus be used, in general, as a simple, readily available convergence accelerator for binding energies. We conclude this section with the following remarks: our procedure for constructing the ABFs gives highly accurate and reliably results for the RI-V approximation. For light elements, one can readily get meV /atom accuracy in total energies and sub-meV /atom accuracy in binding energies, for a wide range of thresholding parameters. For heavy elements, meV /atom accuracy requires tigher thresholding parameters, particulary for the on-site orthonormalization ε orth . However, sub-meV BSSE-corrected binding energy accuracy can always be obtained independent of the choice of thresholding parameters.
Choosing converged yet efficient thresholding parameters thus obviously depends on the element in question. In fact, ε orth can be chosen differently for each element in the same calculation. For RI-V and light elements (Z = 1-10), we employ ε orth = 10 −2 in the following. For heavier elements (Z > 18), we resort to ε orth = 10 −4 , which yields negligibly small errors in total energies even for heavy elements, and ε orth = 10 −3 for elements in between. The additional, system-wide SVD threshold ε svd is set to 10 −4 or tighter for the remainder of this paper. As shown above, its accuracy implications are then negligible as well.

NAO basis convergence for HF, hybrid density functionals, MP2, RPA, and GW
Having established the quality of our ABFs for given orbital basis sets, we next turn to the quality of our actual NAO orbital basis sets for HF, hybrid density functionals, MP2, RPA, and GW calculations. Below we will separate the discussions of self-consistent ground-state calculations (HF and hybrid density functionals) and correlated calculations (MP2, RPA, and GW ). As will be demonstrated below, with our basis prescription, a convergence of the absolute HF total energy to a high accuracy (mev/atom) is possible for light elements. In other cases (heavy elements or correlated methods), the total energy convergence is not achieved at the moment, but we show that energy differences, which are of more physical relevence, can be achieved to a high quality.

HF and hybrid density functional calculations
Here we will demonstrate how well the generic NAO basis set libraries described in section 4.1 and in [64] perform for HF and hybrid density functional calculations. As described earlier, our orbital basis sets contain a functional-dependent minimal basis, composed of the core and valence orbitals of the free atom, and additional functionalindependent optimized basis sets (tiers). Thus, besides the discussion of the convergence behaviour of the generic optimized tiers basis sets, here we will also address the influence of the choice of the minimal basis which is in practice generated by certain atomic solver . For all-electron calculations for molecular systems with a given electronic-structure method, it would be best if the core basis functions were generated from the atomic solver using the same method. In this way, the behaviour of the molecular core wavefunctions in the vicinity of the nuclei would be accurately described at a low price. This is the case for LDA and most GGA calculations in FHI-aims. Similarly, for HF molecular calculations, it would be ideal if the minimal basis was generated by the HF atomic solver. Unfortunately at the moment the HF atomic solver is not yet available in our code, and in practice we resort to the minimal basis generated from other types of atomic solvers. This is not a fundamental problem, and the only price one has to pay is that more additional tiers basis functions are needed to achieve a given level of basis convergence. Nevertheless one should keep in mind that there is a better strategy here and in principle our basis prescription should work even better than what is reported here.
In the following the NAO basis convergence for HF will be examined, and along the way the influence of the minimal basis will be illustrated by comparing those generated by DFT-LDA and Krieger-Li-Iafrate (KLI) [159] atomic solvers. The KLI method solves approximately the exact-exchange optimized-effective-potential (OEP) equation [160,102], by replacing the orbital-dependent denominator in the Green function (at zero frequency) appearing in the OEP equation by an orbital-independent parameter. This in practice reduces the computational efforts considerably without losing much accuracy [161]. The KLI atomic core wavefunctions resemble the HF ones much better than the LDA ones do. As demonstrated below, by moving from the LDA minimal basis to the KLI ones in HF calculations, the abovementioned problem is alleviated to some extent.

Light elements
We first check the convergence behaviour of our NAO basis sets for light elements in HF and hybrid functional calculations. In figure 7 (left panel) the HF total energy of N 2 as a function of increasing basis set size is plotted, starting with two different sets of minimal bases -generated respectively from LDA and KLI atomic solvers. All other basis functions beyond the minimal part (the tiers) are the same for both curves. For comparison, the convergence behaviour of the LDA total energy of N 2 is shown on the right panel for the same basis sets. The horizontal (dotted) lines indicate independently computed GTO reference values using NWChem and the Dunning cc-pV6Z basis set, which gives the best estimate for the HF total energy at the complete basis set limit [162,163].
First, with both types of minimal basis the HF total energy can be systematically converged to within a few meV of the independent GTO reference value. This is reassuring, as we can thus use our standard NAO basis sets in a transferable manner even between functionals that are as different as LDA and HF. Furthermore, it is evident that the KLI-derived minimal basis performs better for HF, and similarly the LDA derived minimal basis performs better for LDA. As mentioned above, this is because the closer the starting atomic core basis functions to the final molecular core orbitals, the faster the overall basis convergence is. If the true HF minimal basis was used, we should expect an even faster basis convergence of the HF total energy, similar to the LDA total-energy convergence behaviour starting with the LDA minimal basis ((blue) circles in the right panel of figure 7). In this case the BSSE in a diatomic molecular calculation should also be vanishingly small since the atomic reference is already accurately converged from the outset with the minimal basis. In practice, the reliance on KLI-derived minimal basis functions instead leads to some small BSSE-type errors in energy differences, as shown below.
In figure 8 we present the NAO basis convergence of the HF binding energy for N 2 as a function of bond distance. Here we start with the KLI minimal basis and then systematically add basis functions from tier 1 to tier 4. Results are shown both without (left) and with (right) a counterpoise correction. A substantial improvement of the binding energy is seen between the tier 1 and the tier 2 basis, with only slight changes beyond tier 2. In the absence of BSSE corrections a slight overbinding is observed for tier 2 and tier 3. As noted above, we attribute this to the fact that in our calculations we used KLI core basis functions as a practical compromise and the atomic reference calculation is not sufficiently converged for the outset. The basis functions from the neighboring atom will then still contribute to the atomic total energy and this leads to non-zero BSSE. The counterpoise correction will cancel this contribution-which is very similar for the free atom and for the molecule-almost exactly. The counterpoise- The two convergence curves correspond to starting minimal basis sets generated by LDA and KLI atomic solvers, respectively. The dotted horizontal line marks the HF and LDA total energy computed using NWChem and the cc-pV6Z basis, which gives a reliable estimate of the basis-set limit within 2 meV [162]. corrected binding energies for tier 2 are in fact almost the same as for tier 4. The latter agrees with results from a GTO cc-pV6Z basis (NWChem) within 1-2 meV (almost indistinguishable in figure 8). Figure 9 demonstrates the same behaviour for a different test case, the binding energy of the water dimer using HF (left) and the PBE0 [30,7] hybrid functional (right). The geometry of the water dimer has been optimized with the PBE functional and a tier 2 basis. For convenience, the binding energy is computed with reference to H 2 O fragments with fixed geometry as in the dimer, not to fully relaxed H 2 O monomers. This is sufficient for the purpose of the basis convergence test here. Detailed geometrical information for water dimer can be found in [64]. In figure 9 the dotted line again marks the NWChem GTO cc-pV6Z reference results. Similar to the case of N 2 , we observe that the HF binding energy is fairly well converged at the tier 2 level, particularly after a counterpoise correction, which gives a binding energy that agrees with the NWChem reference value to within 1-2 meV. The BSSE arising from insufficient core description is reduced for the PBE0 hybrid functional, where only a fraction (1/4) of exact-exchange is included. In practice, HF calculations at the tier 2 level of our NAO basis sets yield accurate results for light elements. Counterpoise corrections help to cancel residual total-energy errors arising from a non-HF minimal basis. However, even without such a correction, the convergence level is already pretty satisfying (the deviation between the black and the red curve in figure 9 is well below 10 meV for tier 2 or higher).

Heavy elements
For heavier elements (Z> 18), the impact of non-HF core basis functions on HF total energies is larger. However, the error again largely cancels in energy differences, as will be shown below. In order to avoid any secondary effects from different scalar-relativistic approximations to the kinetic energy operator, in figure 10 (upper panels) we first compare the convergence of non-relativistic (NREL) HF total energies with NAO basis size for the coinage metal dimers Cu 2 , Ag 2 , and Au 2 at fixed binding distance d=2.5Å. (The experimental binding distances are 2.22Å [164], 2.53Å [165,166] and 2.47Å [164], respectively.) Again, we find that KLI-derived minimal basis HF total energies (upper panels), non-relativistic binding energies (middle panels) and scalar-relativistic (scaled ZORA [158,64]) binding energies (bottom panels) of Cu 2 (left), Ag 2 (middle), and Au 2 (right), at fixed bond length, d=2.5Å. Similar to figure 7, results are shown for two sets of minimal basis generated using LDA and KLI atomic solvers. For clarity the NREL total energies are offset by -89197.12 eV, -282872.42 eV, and -972283.08 eV respectively for Cu 2 , Ag 2 , Au 2 , which correspond to the actual values of the last data points with KLI minimal basis. All binding energies are BSSE-corrected. The dashed horizontal lines in the bottom panels mark the NWChem reference values using aug-cc-pV5Z-PP basis with ECP.
sets are noticeably better converged (lower total energies) than LDA-derived minimal basis sets. In contrast to N 2 , however, absolute convergence of the total energy is here achieved in none of these cases, and the discrepancy increases from Cu (nuclear charge Z= 29) to Au (Z= 79). For comparison, the middle and lower panels of figure 10 show non-relativistic and scalar-relativistic binding energies for all three dimers. The scalar-relativistic treatment employed is the scaled ZORA due to Baerends and coworkers [158] (for details of our own implementation, see [64]). In all three cases, the binding energies are converged to a scale of ∼0.02 eV, at least two orders of magnitude better than total energies. In other words, any residual convergence error due to the choice of minimal basis (LDA or KLI instead of HF) cancel out almost exactly. To compare our prescription to that generally used in the quantum chemstry community where effective core potentials (ECP) are used to describe the core electrons and the relativistic effect, we also marked in figure 10 the reference values computed using NWChem and the aug-cc-pV5Z-PP basis [167,168]. The agreement between our all-electron approach and the GTO-ECP one is pretty decent for Ag 2 and Au 2 , whereas a larger discrepancy of ∼0.03 eV is observed for Cu 2 . For this particular case we suspect the remaining disagreement is an issue with respect to the atomic reference energy for the Cu atom between both codes. More work will be done to fully unravel the point.

MP2, RPA, and GW calculations
In the implementation described here, MP2, RPA, and GW methods require the explicit inclusion of unoccupied single-particle states. As a consequence, noticeably larger basis sets are needed to obtain converged results in these calculations [169,170,171,172,173,174,175,38,86,90]. Much experience has been gained in the quantum chemistry community to construct Gaussian basis sets for correlated calculations [157,176], but for NAOs this is not case. In this section, we show how our standard NAO basis sets perform for MP2, RPA, and GW calculations, for both light and heavy elements. For clarity, we separate the discussions for the convergence of binding energies (in the case of MP2 and RPA) and quasiparticle excitations (in the case of GW and MP2 self-energy calculations). In contrast to the cases of HF and hybrid density functionals, BSSE corrections for RPA and/or MP2 are essential to obtain reliable binding energies. This results directly from the larger basis sets required to converge the MP2 or RPA total energy [169,170,171,172,38,86,90], yielding larger BSSE for finite basis set size. With our standard NAO basis sets, the actual BSSEs in MP2 and RPA (based on the HF reference, denoted as RPA@HF in the following) calculations are plotted in figure 11 for the example of N 2 . The size of BSSEs in these cases is huge and does not diminish even for the pretty large tier 4 basis. It is thus mandatory to correct these errors in MP2 and RPA calculations to get reliable binding energies. As one primary interest in this work is the applicability of standard NAO basis sets for MP2 and RPA, all binding energies presented are therefore counterpoise-corrected. In all HF reference calculations in this session, the KLI minimal basis is used. For RPA, GW , and MP2 self-energy calculations, we use 40 imaginary frequency points on a modified Gauss-Legendre grid (7), which ensures a high accuracy for the systems studied here. figure 12 the BSSE-corrected MP2 and RPA binding energies for N 2 and the water dimer are shown as a function of the NAO basis set size. The dotted line marks reference results computed with FHI-aims and the Dunning aug-cc-pV6Z basis. In the case of MP2, the FHI-aims aug-cc-pV6Z values agree with that of NWChem to within 0.1 meV. Unfortunately, a similar independent reference is not available for RPA, but excellent agreement is also seen with smaller basis sets, for which reference RPA data are available for N 2 [38].  Upon increasing the basis size, the biggest improvement occurs when going from tier 1 to tier 2, with further, smaller improvements from tier 2 to tier 4. For the strongly bonded N 2 the MP2 binding energy at tier 4 level deviates from the aug-cc-pV6Z result by ∼ 120 meV, or ∼ 1% of the total binding energy. For the hydrogen-bonded water dimer, the corresponding values are ∼ 3 meV and ∼ 1.5% respectively. The convergence quality of RPA results with respect to the NAO basis set size is similar. Going beyond our FHI-aims standard NAO basis sets, further improvements arise by adding (ad hoc, as a test only) the diffuse functions from a GTO aug-cc-pV5Z basis set, denoted "a5Z-d" in the following. The results computed using this composite "tier 4 + a5Z-d" basis are shown by the last point in figure 12. The deviation between the tier 4 and aug-cc-pV6Z results is then reduced by more than a factor of two. For the water dimer, for example, "tier 4 + a5Z-d" gives -220.9 meV and -206.5 meV for the MP2 and RPA@HF binding energies, comparable to the quality of the cc-pV6Z basis which yields -221.1 meV and -206.9 meV, respectively. Both then agree with the aug-cc-pV6Z results (-222.3 meV for MP2 and -208.9 meV for RPA@HF) to within ∼ 2 meV. In this context, it is interesting to check if the cut-off radii of our NAO functions have any influence on the convergence behaviour demonstrated above. As described in [64], the NAO basis functions are strictly localized in a finite spatial area around the nuclei, and the extent of this area is controlled by a confining potential. For the default settings used in the above calculations, this potential sets in at a distance of 4Å from the nucleus and reaches infinity at 6Å. The question is what would happen if we reduce or increase the onset radii of this confining potential? The answer to this question is illustrated in figure 13 where basis convergence behaviour for N 2 and the water dimer are shown for three different onset distances of the confining potential. From figure 13 one can see that increasing the onset radius of the confining potential (i.e., enlarging the extent of the NAO basis functions) from the default value (4Å) has little effect on the convergence behaviour for N 2 or (H 2 O) 2 . Upon reducing it, noticeable changes of the results only occur for tier 1 or tier 2 in certain cases, but the overall effect is very small and does not change the general convergence behaviour described above. This finding holds in general for covalent and hydrogen bonds. In practice, the onset radius Results from an independent, Gaussian-type calculation (aug-cc-pV5Z-PP basis set with ECP, NWChem code) are included for comparison. may always be invoked as an explicit convergence parameter-for instance, much more weakly bonded (dispersion bonded) systems benefit from slightly larger radii (5Å -6 A) in our experience. Further details on this can be found in [177].

Binding energies As illustrating examples for light elements, in
We next illustrate the NAO basis convergence for heavy elements, using Au 2 as an example. In figure 14 the MP2 binding energy for the Au 2 dimer as a function of the bond length is plotted for different NAO basis sets. Relativity is again treated at the scaled ZORA level [158,64]. The binding curves shown here demonstrate that the same qualitative convergence behaviour as for our light-element test cases carries over. In essence, significant improvements are gained from tier 1 to tier 2, and basis sets between tier 2 and tier 4 yield essentially converged results. For comparison, we show a completely independent (NWChem calculations) curve with Gaussian "aug-cc-pV5Z-PP" basis sets [167,168]. The resulting binding energy curve yields rather close agreement with our all-electron, NAO basis set results.

Quasiparticle energies
After the discussion of binding energies, we next examine how the GW and MP2 quasiparticle energy levels converge with our NAO basis sets. The G 0 W 0 @HF and MP2 quasiparticle HOMO levels for N 2 and the water dimer are plotted in figure 15 as a function of basis set size. MP2 is denoted here as "MP2-QP" to emphasize that the MP2 self-energy (32) is used, rather than a MP2 totalenergy difference. We again take the results of an "aug-cc-pV6Z" GTO calculation as a reference. The first four data points in each sub-plot correspond to the NAO basis sets. The last point represents the composite "'tier 4 + a5Z-d" basis as described above. Once again, the biggest improvement occurs when going from tier 1 to tier 2. However, the deviation of the HOMO levels from the reference values at tier 2 level is The remaining error is then further reduced by a factor of two by including the diffuse "a5Z-d" part of a 5Z GTO basis set. Accounting for the possible underconvergence of the aug-cc-pV6Z itself (compared to the CBS limit), we expect an overall ∼0.1 eV under-convergence of the composite "'tier 4 + a5Z-d" results given here. This is still an acceptable accuracy, considering the generally known challenge of converging the correlation contribution involving virtual states using local orbital basis functions [169,170,171,172,38,86,90]. Therefore "'tier 4 + a5Z-d" basis sets were used in the benchmark calculations presented in the section 6.

5.2.3.
G 0 W 0 calculations for the benzene molecule G 0 W 0 calculations for molecules have been reported in a number of publications in recent years using various numerical frameworks and computer code packages. In particular in the solid-state community, calculations have been based on plane wave methods together with the supercell approach and pseudopotential approximations. Their advantage is a systematically convergable basis set (plane waves) especially for the unoccupied spectrum. However, the two approximations mentioned above (pseudopotential and supercell) can be drastic [178,179,180,181,182,183,184]. In particular, the Coulomb operator in vacuum is not screened in a standard plane wave approach. As a result, different images of the system interact with one another across supercell boundaries [180,181,182,183,184]. In addition, the slow converence of G 0 W 0 results with the plane-wave cutoff of the unoccupied spectrum is notorious [175,185,186]. Table 3 reports literature G 0 W 0 results for the benzene molecule[187, 128, 188, Table 3. G 0 W 0 HOMO and LUMO values for the benzene (C 6 H 6 ) molecule obtained by FHI-aims and several other numerical approaches as reported in literature. "A.E" and "P.P." in the second column refer to "all electron" and "pseudopotential" respectively.  [186] : plane wave basis augmented with siesta-type localized atomic orbitals c Ref. [128] : plane wave basis extrapolated to infinite energy cutoff d Ref. [188] : plane wave basis plus Lanczos trick to remove the empty states e Ref. [153] f Ref. [187] g Ref. [189] h Ref. [190] : (negative) ionization energy (IE). The vertical IE only differs from this value by 0.01 eV according to the NIST database [191]. i Ref. [192] : (negative) vertical electron affinity (IE) 186,153,189], a particularly often studied case, in comparison to our own results. We focus on G 0 W 0 @LDA and G 0 W 0 @HF for the HOMO and LUMO levels. Based on these results, it is clear that there is a significant degree of scatter, even between results that are ostensibly converged using the same fundamental approximations. Our own results for NAO "tier 3", "tier 4", and "tier 4 + a5Z-d" basis sets suggest internal convergence at the "tier 4" level: −9.06 eV for the HOMO, and 0.94 eV for the LUMO in G 0 W 0 @LDA, compared to −9.64 eV and 1.51 eV in G 0 W 0 @HF. The LUMO values are unbound and can be interpreted as experimental resonance energies (here taken from the tabulated negative vertical electron affinity). In either case (HOMO or LUMO, G 0 W 0 @LDA or G 0 W 0 @HF) the results are not far from the experimental values.
The same cannot be said for the comparison between the different numerical implementations. For instance, the HOMO values in G 0 W 0 @LDA range from −8.78 eV (small numeric basis set and pseudopotentials) to −9.88 eV on a real-space grid. Even within the plane wave based approaches, the values range from −9.03 eV to −9.40 eV. Similar discrepancies arise for G 0 W 0 @HF, and for the LUMO values. Within the scatter evidenced by Table 3, we believe that our own results possess some merit, as the system (i) is isolated (no supercell), (ii) is treated fully all-electron, and (iii) any residual convergence issue of the unoccupied spectrum should be exposed by the diffuse Gaussian basis functions ("a5Z-d").

Basis set converged benchmark data for the G2 and S22 molecular test sets
Our final section utilizes the preceding methodologies and techniques to provide accurate all-electron results for molecular test sets close to basis set convergence. We cover vertical ionization energies and binding energies at different levels of theory, for well-defined, published molecular geometries.
6.1. Vertical ionization energies from HF, MP2, and G 0 W 0 methods The G 0 W 0 method has been used to calculate the single-particle properties for various small-and medium-size molecules with considerable success [193,194,195,187,196,197,198,199,200,128,201,202,189,153]. The dependence of G 0 W 0 on the starting point has also been looked at in the past [60,61,197,202]. Here, using a selected collection of atoms and molecules from the G2 ion test set for ionization energies [79], we aim to establish the overall performance of G 0 W 0 for computing IEs close to the basis set limit, and systematically examine its dependence on the starting point. The selection of molecules is based on the availability of experimental geometries and experimental vertial IEs. With vertical IEs we denote ionization energies at fixed geomerty, i.e. no structural relaxation after excitation. This quantity is directly comparable to the GW and MP2-QP and quasiparticle energies. We will also assess the MP2-QP approach for determining IEs since this method was used in quantum chemistry in past decades, but direct comparisons in the literature are scarce [44]. Finally, results for IEs determined by MP2 total energy differences (denoted here simply as "MP2" in contrast to "MP2-QP") between the neutral atoms/molecules and the corresponding positively charged ions are also presented. We note that IEs given by MP2-QP essentially correspond to the MP2 total energy difference if the HF orbitals of the neutral systems were used in the ionic calculation (for a discussion, see, e.g., [44]). All calculations were carried out at the experimental geometries. As mentioned above, the composite "tier 4 + a5Z-d" basis set is used for most of the elements except for a few cases (e.g., rare gas atoms) where the tier 4 NAO basis is not available. In these cases the full aug-cc-pV5Z functions are added to the NAO tier 2 basis. On average we expect the chosen basis setup to guarantee a convergence of the HOMO quasiparticle levels to ∼0.1 eV or better for the calculations presented in the following.
In figure 16 we present histograms of the error distributions given by HF, MP2-QP, G 0 W 0 @HF and G 0 W 0 @PBE0 for this database. The actual IE values are presented in 7. Compared to HF, one can see that the deviations from the experimental values are much smaller in G 0 W 0 and MP2-QP. On average HF tends to overestimate IEs. This trend is corrected by MP2 and G 0 W 0 . Using the same HF reference, the magnitude of   Figure 16. Histograms of the error distribution for IEs for a set of 50 atoms and molecules calculated with HF, G 0 W 0 @HF, G 0 W 0 @PBE, G 0 W 0 @PBE0, MP2-QP, and MP2.
the correction is smaller in G 0 W 0 than in MP2-QP, due to the renormalization effect coming from the screened Coulomb interaction in the GW self-energy. Concerning the starting-point dependence of the G 0 W 0 method, G 0 W 0 based on HF gives too large IEs on average, whereas G 0 W 0 based on PBE does the opposite. Consequently, G 0 W 0 based on the hybrid density functional PBE0 appears to be the best compromise, although a slight underestimation of the IEs is still visible. Furthermore, comparing MP2-QP with MP2 shows that the two approaches yields very similar results, implying that for the light elements reported here orbital relaxation effects are not significant [203]. Table 4 summarizes the error statistics. For the subset of atoms and molecules MP2-QP gives the smallest mean error (ME) and G 0 W 0 @PBE0 the smallest mean absolute error (MAE).

Benchmark MP2 and RPA results for the G2-I atomization energies
The 55 atomization energies from the original G2-I set (Table III of [204]) serve as a good benchmark database for ground-state total-energy methods for covalently bonded systems. In Table 5 we present our MP2 and RPA atomization energies for the G2-I set as computed using FHI-aims and the "tier 4 + a5Z-d" basis set. The geometries used in the calculations are those determined by Curtiss et al. [204], i.e., all-electron MP2 optimization with a Gaussian 6-31G* basis set. The reference data, taken from [205], is derived from experiment and corrected for zero-point vibrations.
Our MP2 and RPA results are in good agreement with those reported by Feller and Peterson [206] and Paier et al. [90], respectively, both obtained with GTO basis sets and extrapolated to the complete basis set limit. Table 5 demonstrates that MP2 performs best for these small covalently bonded molecules. The RPA approach, on the other hand, which has a broader applicability (e.g. bond breaking and/or metals where MP2 fails), exhibits significant underbinding. For RPA we also observe that KS-PBE provides a much better reference than HF. Recently computational approaches to overcome the underbinding behaviour of the standard RPA scheme have been proposed. These comprise contributions from single excitations [40] and second-order screened exchange [110,90].

Benchmark MP2 and RPA binding energies for the S22 molecular set
The S22 molecular set proposed by Jurečka et al. [80] has become a standard benchmark database for testing the accuracy of existing and newly developed methods for the description of weak interactions. S22 represents an "unbiased" set in the sense that it contains molecules of different bonding nature (7 hydrogen bonded, 8 dispersion bonded, and 7 with mixed bonding character) and of different size (ranging from small ones like the water dimer to relatively large ones like the Adenine-thymine dimer containing 30 atoms). The MP2 and CCSD(T) interaction energies for this set of molecules were already computed by Jurečka et al. and extrapolated to the Gaussian complete basis set (CBS) limit. A more consistent and accurate extrapolation for the CCSD(T) values was recently carried out by Takatani et al. [207], which we therefore adopt as reference here. While MP2 calculations for S22 are common, RPA benchmark calculations for the whole set have not been performed, yet. RPA-type calculations for S22 have recently been reported [98,40,208]. The agreement between the data from different groups is not perfect, and the basis incompleteness could be an issue. Of our own RPA numbers only the MAEs and MAPEs have been reported previously [40]. Now in Table 6 the actual MP2, RPA@HF, and RPA@PBE binding energies for these molecules as computed using FHI-aims and the composite "tier 4 + a5Z-d" basis set. With our basis setup, the MP2 binding energies are underestimated by ∼4.5 meV (2% on a relative scale) compared to the MP2/CBS results reported in [80]. We expect a similar convergence of the RPA numbers based on our basis set convergence tests shown in previous sections. Compared to the CCSD(T) reference data, RPA@PBE, which nowadays dominates  RPA-type production calculations, systematically underestimates all of the three weak bonding categories and gives a MAE of 39 meV and a MAPE of 16%. If instead HF is used as a starting point (i.e., as for MP2), the description of hydrogen bonding improves, while the description of dispersion bonding worsens. The overall MAE for RPA@HF is 51 meV and the MAPE 25%. We are thus faced with the conundrum that RPA can describe the weak interactions that are beyond the reach of LDA, GGA, and hybrid functionals, but the accuracy of the two standard RPA schemes is not spectacular. We have analyzed the origin of this underbinding behaviour in [40] and proposed a simple solution to overcome this problem.

Conclusions and outlook
To summarize, we have presented a resolution of identity framework for the two-electron Coulomb operator that allows efficient, accurate electronic structure computations based on HF, hybrid functionals, MP2, RPA, and GW based on the flexible basis function form of NAOs. We have shown that NAO basis sets as implemented in FHI-aims are a competitive choice for approaches involving exact-exchange and/or non-local correlation terms, with rather compact basis sets sufficing for essentially converged results. Our simple "on-the-fly" scheme to construct the auxiliary basis functions using Gram-Schmidt orthonormalization of the "on-site" products of single-electron atomic orbitals gives a natural, accurate representation of the two-electron Coulomb operator for practical calculations. Taken together, our framework paves the way for an extended usage of NAOs in more advanced computational approaches that go beyond LDA and GGAs. Specifically, we have applied the G 0 W 0 and MP2 quasiparticle approaches to compute the vertical IEs of a set of small molecules, and the RPA method to compute Table 6. Binding energies (in meV) of the S22 molecular set [80] calculated with MP2, RPA@HF, and CCSD(T). MP2 and RPA results have been obtained with FHI-aims and "tier 4 + a5Z-d" basis sets. The reference CCSD(T) results are from [207].
particular towards the limit of very large systems: the benefit of compact basis set size at a given accuracy level should help tackle problem sizes that, otherwise, could not be done.
technology that is used in many quantum-chemical applications for the exchangecorrelation matrix in DFT [154,65]. The integration grid points r = r(a, s, t) are uniquely specified by the atomic center a, the radial shell number s, and angular point t. w(s, t) is the corresponding integration weight. Details of our own implementation (FHI-aims) are given in [64,155]. Since these are true three-center integrals, we restrict the integration domain for a particular integral element to the grids associated with the atoms on which the basis functions in question are centered. For instance, denoting the respective atoms by a 1 , a 2 , a 3 , the three-center integrals in (66) can then be discretized as (ij|µ) = a=a 1 ,a 2 ,a 3 s,t p 3 (a, r)w(s, t)φ i (r)φ j (r)Q µ (r), (. where p 3 (a, r) is a three-center partition function that satisfies a=a 1 ,a 2 ,a 3 p 3 (a, r) = 1 everywhere in the overlapping region of the three functions, and is zero otherwise. The underlying numerical grids can in principle be increased up to arbitrary accuracy if needed. The two-center integrals (67) can be performed in a similar fashion using overlapping grids, or with the spherical Bessel transform techniques explained below.
The three-dimensional integral (.6) can be separated by expanding the plane wave e ik·R in spherical harmonics and spherical Bessel functions The separation yields f (r)g(r − R) dr = (2π) in (.12) can be calculated efficiently using recursion formulae [213]. They are nonzero only for L = |l −l ′ |, |l −l ′ | + 2, . . . , (l + l ′ ). If two atom-centered functions do not overlap, the overlap integrals I L (R) of course vanish. For any given distance R, the radial integrals I L (R) in (.11) can be calculated directly using the trapezoidal rule on the logarithmic grid when j L (kR) is evaluated in logarithmic Fourier space as described in the next section. If integrals of the same atom-centered functions for many differing distances are needed, one can compute these more efficiently by interpreting (.11) as an SBT ofP (k) =f (k)g(k). By applying the logSBT, and interpolate for all needed distances R, one can obtain all the integrals at tight-binding cost, meaning that efficient recursion formula (together with spline evaluations) can be used intead of evaluating each integral numerically.
Coulomb interactions of atomic functions can be calculated with comparable ease where the integrand in (.11) is multiplied with the Coulomb kernel 4π/k 2 V L (R) = 4π 2 π ∞ 0 j L (kR)f (k)g(k) k 2 k 2 dk. (.14) The Coulomb interaction generally does not vanish even if the two charge densities do not overlap. However, it has a simple multipolar behaviour and explicit integration of (.14) can thus be avoided. The function V L (R) can formally be interpreted as the far field of a charge distribution of angular momentum L whose radial part P (r) is given by its SBTP (k) =f (k)g(k). Therefore, it only depends on the multipole moment of P , which is encoded in its limiting behaviour for small k. From this, it can be shown that V L (R) vanishes for all L < l + l ′ . For L = l + l ′ we get with (2n − 1)!! = 1 · 3 · · · (2n − 1) and with multipole moments (.16) Therefore, the Coulomb interaction of non-overlapping functions depends only on the product of their multipole moments and can also be calculated at tight-binding cost.
As a side remark we emphasize that two ABFs do not interact via the Coulomb interaction if they do not overlap and at least one of the two multipole moments is zero. As shown by Betzinger in [217], all but one of the ABFs for a given atom and a given angular momentum lm can then be chosen to be multipole free by means of a suitable unitary transformation.

Appendix A.4. Logarithmic spherical Bessel transform
This section describes our implementation of the logSBT algorithm. For a more extensive description we refer the interested reader to Talman [216] and Hamilton [215].
The SBT as defined in (.8) can be written as the integral over a kernel J (kr) and a right-hand-side F (r): .

(.17)
The choice of the power bias parameter 0 ≤ α ≤ 3 is crucial for the numerical accuracy and stability of this method and will be discussed further below. The basic idea of the logSBT is that in logarithmic coordinates (ρ = log r and κ = log k) the kernel reads J (κ + ρ) and (.17) turns into a convolution, which can be efficiently calculated using fast Fourier transorms (FFTs). Please note that dr/r = dρ in (.17). As pointed out by Hamilton [215], this procedure is exact if both F (ρ) = e (3−α)ρ f (e ρ ) and the corresponding SBT termF (κ) = e ακf (e κ ) are periodic in logarithmic space and analytic expressions for the logarithmic Fourier transform of the kernel are used. Periodicity can be achieved, e. g., by choosing α near 1.5 and using a sufficiently wide logarithmic grid. Under these circumstances, both F (ρ) andF(κ) smoothely drop to zero on both ends, which therefore can safely be connected.
Unfortunately, the scaling factor k −α turns out to be quite problematic. By design of the algorithm, the absolute error ofF (κ) = k α f (k) before the final scaling is typically of the order of machine precision, i. e., about 10 −15 . This is true even if the magnitude of the exact value is much smaller than that. After scaling, however, the absolute error can get arbitrarily large because k can be very small on a wide logarithmic grid.
Talman [216] circumvents this problem by using two separate α for small and large k and joining the two results where they differ least. The small-α (α = 0) calculation cannot be done assuming periodicity for l = 0 because F (ρ) does not decay to zero for ρ → −∞. Therefore, a trapezoidal rule is used for the integration, which works well for small k where the Bessel function is smooth on a logarithmic scale. This does not break with the spirit of logSBT, because the trapezoidal rule can be formulated using FFTs, too.
In order to avoid the second transform, we take a different approach. In practice, one needs the SBT only for one kind of integral and it is sufficient to calculate k αf (k) to high absolute accuracy for a single given α, which can be used as power bias for the transform.
This works well for all cases but α = 0 and l = 0. Here, we cannot simply resort to the trapezoidal rule because it is invalid for high k where the Bessel function oscillates rapidly. Instead, we separate j 0 (e τ ) into a smooth part proportional to erfc(τ /∆τ 0 ) and a properly decaying rest. We use the first part for a trapezoidal rule and the second part for the log-periodic algorithm. Fortunatly, these two schemes differ only in the way the kernel is constructed so that the two kernels can simply be added up to a "hybrid" kernel. The actual transform is not affected and numerically not more expensive than an ordinary logSBT.
Just like Talman [216] we double the domain during the transforms for l = 0 and l = 1 in order to avoid the need for large domains for a proper decay behaviour.

Appendix B. Ionization energies of the a set of atoms and molecules
In Table 1 the individual numbers for the vertical ionization potentials for 50 atoms and molecules as computed with 6 different computational approaches are presented. The calculations are performed with FHI-aims and "tier +a5Z-d" basis set.

Appendix C. Modified Gauss-Legendre grid
For the integrals over the imaginary frequency axis (e.g., for the RPA correlation energy, equation (60) ), we use a modified Gauss-Legendre quadrature. The Gauss-Legendre quadrature provides a way to numerically evaluate an integral on the interval [−1 : 1] where x i and w i are the integration points and the corresponding weights, respectively. For our purposes a transformation procedure is applied to map the integration range from [−1 : 1] to [0 : ∞] whereby the x i and w i have to be changed accordingly. Table 1. The vertical ionization potentials (in eV) for 50 atoms and molecules (taken from G2 ion test set [79]) calculated with HF, MP2, MP2-QP, G 0 W 0 @HF, G 0 W 0 @PBE0, and G 0 W 0 @PBE in comparison to the experimental values, taken from the the NIST database [191]. The mean absolute errors (MAE) for the three approaches are also shown.
with x 0 set to 0.5. The weights for the tranformed grid are then given bỹ This modified Gauss-Legendre scheme allows a quick convergence of the frequency integration with a relatively small number of frequency points. In our implementation, a 40-point grid gives micro-Hartree total energy accuracy for the systems investigated in this work.