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

HF theory, hybrid density functionals and MBPT (MP2, RPA and GW) are all approximate ways to solve the interacting many-electron Hamiltonian

$$\begin{equation} \skew3\hat{H} = \sum_{i=1}^{N_{\rm e}}\left[-\frac{1}{2} \nabla^2_i + v_i^{\rm ext} \right] + \sum_{i < j}^{N_{{\rm e}}} v^{\rm ee}_{ij} , \end{equation} \tag{ 1 }$$

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 (the 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 $\skew3\hat{H}^0$ that yields an approximate ground-state wave function |Φ₀〉:

$$\begin{eqnarray} \hphantom{|\Phi_0 \rangle}\skew3\hat{H}^{0} &=& \sum_{i=1}^{N_{{\rm e}}} \hat{h}^{0}_i = \sum_{i=1}^{N_{{\rm e}}} \left[ -\frac{1}{2} \nabla^2_i + v_i^{\rm ext} + v^{\rm MF}_i \right], \nonumber\\ \skew3\hat{H}^{0} |\Phi_0 \rangle & =& E_0^{(0)}|\Phi_0 \rangle. \end{eqnarray} \tag{ 2 }$$

The '(0)' in E⁽⁰⁾₀ implies the fact that this is the ground-state energy of the mean-field Hamiltonian. A suitable $\skew3\hat{H}^{0}$ should be solvable with relative ease. |Φ₀〉 is a single Slater determinant formed from the lowest N_e single-particle spin-orbitals determined by

$$\begin{equation} \hat{h}^0|\psi_{n\sigma}\rangle = \epsilon_{n\sigma}|\psi_{n\sigma}\rangle, \, \end{equation} \tag{ 3 }$$

where σ denotes the spin index and $\hat {h}^0$ is the effective single-particle Hamiltonian noted in the brackets of (2). The form of (2)–(3) is, of course, precisely that of the 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

• $\skew3\hat{H}$ (the many-electron Hamiltonian),
• the mean-field Hamiltonian $\skew3\hat{H}^{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 $\skew3\hat{H}$),
• the effective single-particle Hamiltonian $\hat {h}^0$, which generates the single-particle orbitals ψ_nσ and orbital energies _nσ.

In MBPT, one starts from $\hat {h}^0$ and the associated eigenenergies and eigenfunctions to systematically approximate the properties of $\skew3\hat{H}$, e.g. its true ground-state energy E₀ in MP2 or RPA or its single-particle excitations in GW. The interacting many-electron Hamiltonian $\skew3\hat{H}$ is partitioned into a mean-field Hamiltonian $\skew3\hat{H}^{0}$ as given by (2) and an interaction term $\skew3\hat{H}'$,

$$\begin{eqnarray} \skew3\hat{H} &= & \skew3\hat{H}^0 + \skew3\hat{H}' \nonumber \\ \skew3\hat{H}' &=& \sum_{i<j}^{N_{{\rm e}}}\frac{1}{| {\bf{r}}_i- {\bf{r}}_j|} - \sum_{i=1}^{N_{{\rm e}}} v^{\rm MF}_i. \end{eqnarray} \tag{ 4 }$$

In the remainder of this section we collect the basic formulae that define the mean-field 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.

In HF theory the ground-state wave function of the Hamiltonian in (1) is approximated by a single Slater determinant |Φ₀〉 and E⁽⁰⁾₀ is obtained by a variational optimization, leading to

$$\begin{equation} \fl \langle {\bf r}|\hat{f}|\psi_{n\sigma}\rangle= \left[-\frac{1}{2} \nabla^2 + v_{\rm ext}({\bf{r}}) + v^{\rm h}({\bf{r}}) \right] \psi_{n\sigma}({\bf{r}})+ \int {\rm d} {{\bf{r}}'} \Sigma^{\rm x}_{\sigma}({\bf{r}}, {{\bf{r}}'}) \psi_{n\sigma} ({{\bf{r}}'}) = \epsilon_{n\sigma} \psi_{n\sigma}({\bf{r}}). \end{equation} \tag{ 5 }$$

Here $\skew3\hat{f}$ denotes the HF single-particle Hamiltonian, and v^h(r) is the Hartree potential,

$$\begin{equation} v^{\rm h}({\bf{r}}) = \int \frac{n({{\bf{r}}'})}{| {\bf{r}}- {{\bf{r}}'}|} {\rm d} {{\bf{r}}'} \, \end{equation} \tag{ 6 }$$

with the electron density

$$\begin{equation} n({\bf{r}}) = \sum_{n\sigma}^{\rm occ} |\psi_{n\sigma}({\bf{r}})|^2 , \end{equation} \tag{ 7 }$$

and Σ^x_σ is the non-local, exact-exchange potential

$$\begin{equation} \Sigma^{\rm x}_{\sigma}({\bf{r}}, {{\bf{r}}'}) = - \sum_n^{\rm occ} \frac{\psi_{n\sigma} ({\bf{r}})\psi_{n\sigma}^{\ast}({{\bf{r}}'})}{| {\bf{r}}- {{\bf{r}}'}|} . \end{equation} \tag{ 8 }$$

Equations (5)–(8) form a set of nonlinear 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 mean-field potential v^MF_i in (2) and (4). The HF wavefunction |Φ₀〉 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

$$\begin{eqnarray*} &&E_{\rm HF} = \langle \Phi_0 |\skew3\hat{H}^{0} + \skew3\hat{H}'| \Phi_0 \rangle = \sum_{n\sigma}^{\rm occ} \epsilon_{n\sigma} - E^{\rm h} - E^{\rm x}, \end{eqnarray*}$$

where the Hartree energy E^h and exact-exchange energy E^x are given, respectively, by

$$\begin{eqnarray} E^{\rm h} &=& \frac{1}{2} \int {\rm d} {\bf{r}} \: n({\bf{r}}) v_{\rm h}({\bf{r}}), \nonumber \\ E^{\rm x} &=& \frac{1}{2} \sum_{n\sigma}^{\rm occ} \int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \psi_{n\sigma}^{\ast}({\bf{r}}) \Sigma^{\rm x}_{\sigma}({\bf{r}}, {{\bf{r}}'}) \psi_{n\sigma}({{\bf{r}}'}) \nonumber \\ &=& -\frac{1}{2} \sum_{mn\sigma}^{\rm occ} \int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \frac{\psi_{n\sigma}^{\ast}({\bf{r}})\psi_{m\sigma}({\bf{r}}) \psi_{m\sigma}^{\ast}({{\bf{r}}'})\psi_{n\sigma}({{\bf{r}}'})}{| {\bf{r}}- {{\bf{r}}'}|} . \end{eqnarray} \tag{ 9 }$$

In general, the single-particle spin-orbitals ψ_nσ(r) are expanded in terms of a set of basis functions {φ_i(r)}

$$\begin{equation} \psi_{n\sigma}({\bf{r}}) = \sum_{i} c^{i}_{n\sigma} \varphi_i({\bf{r}}), \end{equation} \tag{ 10 }$$

where cⁱ_nσ are the expansion coefficients. In terms of these basis functions, the HF effective potential can be expressed in a matrix form

$$\begin{eqnarray} V^{\rm HF}_{ij,\sigma} &=& \int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \varphi_i({\bf{r}}) \left[ v^{\rm h}({\bf{r}})\delta({\bf{r}}- {{\bf{r}}'})+\Sigma^{\rm x}_\sigma({\bf{r}}, {{\bf{r}}'})\right] \varphi_j({{\bf{r}}'}) \nonumber \\ &=& v^{\rm h}_{ij,\sigma}+\Sigma^{\rm x}_{ij,\sigma}, \end{eqnarray} \tag{ 11 }$$

where in particular the exact-exchange matrix is given by

$$\begin{equation} \Sigma^{\rm x}_{ij,\sigma} = \sum_{k l} (ik|lj)D_{kl,\sigma}. \end{equation} \tag{ 12 }$$

In (12) D_kl,σ is the density matrix

$$\begin{equation} D_{kl,\sigma} = \sum_{n}^{\rm occ} c_{n\sigma}^{k}c_{n\sigma}^{l\ast} , \end{equation} \tag{ 13 }$$

and (ij|kl) is the shorthand notation of quantum chemistry for four-center two-electron integrals

$$\begin{equation} (ij|kl) = \int\!\!\!\int \frac{\varphi_i({\bf{r}})\varphi_j({\bf{r}})\varphi_k({{\bf{r}}'})\varphi_l({{\bf{r}}'})} {| {\bf{r}}- {{\bf{r}}'}|} {\rm d} {\bf{r}}\, {\rm d} {{\bf{r}}'}. \end{equation} \tag{ 14 }$$

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₀ and ground-state density n₀(r). In practice, the XC energy functional and potential have to be approximated. The effective single-particle orbitals either from HF or from approximate KS-DFT are convenient starting points for MBPT.

Assuming that $\skew3\hat{H}'$ is benign and can be treated as a perturbation on top of $\skew3\hat{H}^0$, the ground-state energy for the interacting system can be obtained using the Rayleigh–Schrödinger perturbation theory (RSPT) or more precisely the 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. This guarantees the size-extensivity of the theory, i.e. the total energy scales correctly with the system size. The Møller–Plesset (MP) perturbation theory is a special case of RSPT [62], where the reference Hamiltonian $\skew3\hat{H}^0$ is the HF Hamiltonian $\skew3\hat{H}_{\mathrm { HF}}=\sum _i \skew3\hat{f}_i$. Terminating the expansion at second order gives the MP2 theory, in which the (second-order) correlation energy is given by

$$\begin{equation} E^{(2)}_0 = \sum_{k > 0} \frac{|\langle \Phi_k|\skew3\hat{H}'|\Phi_0 \rangle|^2}{E^{\rm (0)}_{0} - E^{\rm (0)}_{k}}. \end{equation} \tag{ 15 }$$

Here Φ_k are the Slater determinants representing the excited states of $\skew3\hat{H}^0=\skew3\hat{H}^{HF}$, and E⁽⁰⁾_k are the corresponding excited-state energies. $\skew3\hat{H}'$ is given by (4) with v^MF = v^HF. Among all possible excited-state configurations Φ_k, only double excitations contribute in (15). This is because singly excited Φ_k do not couple to the ground state Φ₀ (Brillouin's theorem [44] for the HF reference), whereas even-higher-excited configurations (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,

$$\begin{equation} \fl E^{(2)}_0 = \frac{1}{2} \sum_{mn}^{\rm occ}\sum_{ab}^{\rm unocc}\sum_{\sigma,\sigma^\prime} (ma,\sigma|nb,\sigma^\prime) \left[ \frac{ (am,\sigma|bn,\sigma^\prime)- (bm,\sigma|an,\sigma^\prime)\delta_{\sigma\sigma^\prime}} {\epsilon_{m\sigma} + \epsilon_{n\sigma^\prime} -\epsilon_{a\sigma}-\epsilon_{b\sigma^\prime}} \right], \end{equation} \tag{ 16 }$$

where (ma,σ|nb,σ') are two-electron Coulomb repulsion integrals for molecular orbitals

$$\begin{equation} (ma,\sigma|nb,\sigma^\prime)=\int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \frac{\psi_{m\sigma}^{*}({\bf{r}})\psi_{a\sigma}({\bf{r}}) \psi_{n\sigma^\prime}^{*}({{\bf{r}}'})\psi_{b\sigma^\prime}({{\bf{r}}'})}{| {\bf{r}}- {{\bf{r}}'}|} . \end{equation} \tag{ 17 }$$

The two terms in (16) correspond to the second-order Coulomb (direct) and second-order exchange energy, respectively.

MP2 corresponds to the second-order term in a perturbation theory where the perturbation expansion is essentially based on the bare Coulomb operator (with the HF effective potential subtracted). As such, the MP2 correlation energy diverges for the 3-dimensional 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, 39, 85–98] is the RPA [16, 35, 36, 53, 99] that in the context of MBPT corresponds to an infinite summation of 'ring' diagrams. The choice of the non-interacting reference Hamiltonian $\skew3\hat{H}^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 [16, 53, 99] 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 OEP approach [100–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 χ⁰ on the imaginary frequency axis

$$\begin{eqnarray} E_{\rm c}^{\rm RPA} &=& \frac{1}{2\pi}\int_0^\infty {\rm d}\omega \,{\rm Tr} \left[{\rm ln} \left(\varepsilon({\rm i}\omega)\right) + \left(1-\varepsilon({\rm i}\omega)\right) \right] \nonumber \\ &=& \frac{1}{2\pi}\int_0^\infty {\rm d}\omega \,{\rm Tr} \left[{\rm ln}\left(1-\chi^0 ({\rm i}\omega)v\right) + \chi^0 ({\rm i}\omega)v\right] \nonumber \\ &=& - \frac{1}{2\pi}\int_0^\infty {\rm d}\omega \sum_{n=2}^{\infty}\frac{1}{n}{\rm Tr} \left[(\chi^0 ({\rm i}\omega)v)^n \right]. \end{eqnarray}\noindent \tag{ 18 }$$

The real-space (Adler–Wiser [104, 105]) representation of χ⁰ reads

$$\begin{eqnarray} \chi^0({\bf{r}}, {{\bf{r}}'},{\rm i}\omega)&=&\langle {\bf{r}} |\chi^0({\rm i}\omega) | {{\bf{r}}'} \rangle \nonumber \\ & =& \sum_{\sigma}\sum_m^{\rm occ}\sum_a^{\rm unocc} \frac{\psi_{m\sigma}^{\ast}({\bf{r}})\psi_{a\sigma}({\bf{r}})\psi_{a\sigma}^{\ast}({{\bf{r}}'}) \psi_{m\sigma}({{\bf{r}}'})} {{\rm i}\omega - \epsilon_{a\sigma} + \epsilon_{m\sigma} } + {\rm c.c.}, \end{eqnarray}\noindent \tag{ 19 }$$

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 χ⁰ through

$$\begin{equation} \varepsilon({\bf{r}}, {{\bf{r}}'}, {\rm i}\omega) = \delta({\bf{r}}- {{\bf{r}}'}) - \int {\rm d} {{\bf{r}}^{\prime\prime}} v({\bf{r}}, {{\bf{r}}^{\prime\prime}})\chi^{\rm 0}({{\bf{r}}^{\prime\prime}}, {{\bf{r}}'},{\rm i}\omega) . \end{equation}\noindent \tag{ 20 }$$

Using (17) and (19), the lowest-order term in (18) can be expressed as

$$\begin{equation} \fl {\rm Tr}\left[\chi^0({\rm i}\omega)v\right]=\int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \chi^0({\bf{r}}, {{\bf{r}}'},{\rm i}\omega) v({{\bf{r}}'}, {\bf{r}}) = \sum_\sigma \sum_m^{\rm occ}\sum_a^{\rm unocc} \frac{(ma,\sigma|am,\sigma)} {{\rm i}\omega-\epsilon_{a\sigma}+\epsilon_{m\sigma}} + {\rm c.c.} \end{equation}\noindent \tag{ 21 }$$

Higher-order terms in (18) follow analogously. There is thus a straightforward path to compute χ⁰ 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 as (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–40, 85–89, 91–98, 106–109]. For instance, when the so-called single-excitation and second-order screened exchange contributions [90, 110] are added to EX + cRPA, the resulting accuracy is impressive [40, 41].

Besides ground-state energies, one is often interested in the properties of electronically excited states. Part of this information is, in principle, accessible by taking the difference in the total energies of the N-electron system and the N ± 1-electron systems using approaches discussed above. In practice, this approach mainly works for computing core-level excitations and/or the first IE and electron affinity for (small) finite systems. In contrast, Green function techniques are more convenient and powerful in 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

$$\begin{equation} G({\bf{r}},t; {{\bf{r}}'}, {t^\prime}) = -{\rm i}\langle N |\hat{T} \hat{\psi}({\bf{r}},t)\hat{\psi}^{\dagger}({{\bf{r}}'}, {t^\prime}) |N \rangle, \end{equation} \tag{ 22 }$$

where $|N\rangle =|\Psi _0({\bf{r}}_1,\dots , {\bf{r}}_N)\rangle$ denotes the interacting ground-state wave function of an N-electron system (solution to (1)). $\hat {\psi }({\bf{r}},t)$ and $\hat {\psi }^{\dagger }({{\bf{r}}'}, {t^\prime })$ are field operators in the Heisenberg picture that annihilate and create an electron at space–time point (r,t) and (r',t'), respectively. $\hat {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

$$\begin{equation} \fl \left[\omega +\frac{1}{2}\nabla^2 -v_{\rm ext}({\bf{r}}) - v^{\rm h}({\bf{r}}) \right] G({\bf{r}}, {{\bf{r}}'};\omega) - \int {\rm d} {{\bf{r}}^{\prime\prime}} \Sigma({\bf{r}}, {{\bf{r}}^{\prime\prime}},\omega)G({{\bf{r}}^{\prime\prime}}, {{\bf{r}}'};\omega) = \delta({\bf{r}}- {{\bf{r}}'}). \end{equation}\noindent \tag{ 23 }$$

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 the diagrammatic Feynman–Dyson perturbation theory. The perturbation series is built on a non-interacting Green function G⁰(r,r',ω) that corresponds to the non-interacting single-particle Hamiltonian $\hat {h}^0$. With single-particle orbitals ψ_nσ(r) and orbital energies _nσ of $\hat {h}^0$, one has

$$\begin{equation} G^{\rm 0}_{\sigma}({\bf{r}}, {{\bf{r}}'},\omega)=\sum_n \frac{\psi_{n\sigma}({\bf{r}})\psi^\ast_{n\sigma}({{\bf{r}}'})}{\omega-\epsilon_{n\sigma}-{\rm i}\eta~{\rm sgn}(\epsilon_{\rm F}-\epsilon_{n\sigma})}, \end{equation}\noindent \tag{ 24 }$$

where _F is the Fermi energy and η a positive infinitesimal. For a given G⁰, corrections to the single-particle excitation energies can be computed from approximate perturbative expansions of Σ(r,r';ω). Examples are the GW method and the second-order approximations discussed below.

In the GW approximation proposed by Hedin [54], the self-energy assumes the form

$$\begin{equation} \Sigma^{{\rm GW}}_\sigma({\bf{r}}, {{\bf{r}}'},\omega) = \frac{{\rm i}}{2\pi}\int {\rm d} \omega G_\sigma({\bf{r}}, {{\bf{r}}'}, \omega+\omega^\prime) W({\bf{r}}, {{\bf{r}}'}, \omega^\prime)\,{\rm e}^{{\rm i}\omega\eta}. \end{equation}\noindent \tag{ 25 }$$

Here W is the screened Coulomb potential at the RPA level

$$\begin{equation} W({\bf{r}}, {{\bf{r}}'}, \omega) = \int {\rm d} {{\bf{r}}^{\prime\prime}} \varepsilon^{-1}({\bf{r}}, {{\bf{r}}^{\prime\prime}}, \omega) v({{\bf{r}}^{\prime\prime}}, {\bf{r}}) , \end{equation}\noindent \tag{ 26 }$$

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⁰W⁰) based on a fixed, non-interacting reference state (DFT with popular functionals or HF) are often performed. With χ⁰ computed from the non-interacting Green function in (24), the ensuing W can be expanded in powers of χ⁰v:

$$\begin{equation} W=v+v\chi^0 v + v\chi^0 v\chi^0 v + \cdots. \end{equation}\noindent \tag{ 27 }$$

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⁰W⁰ self-energy is obtained from (25), the corrections to the single-particle orbital energies are given by

$$\begin{equation} \epsilon_{n\sigma}^{G^0W^0} = \epsilon_{n\sigma} + \left\langle \psi_{n\sigma}| \Sigma_\sigma^{G^0W^0} \left(\epsilon_{n\sigma}^{G^0W^0}\right) - v^{\rm xc}| \psi_{n\sigma}\right\rangle , \end{equation}\noindent \tag{ 28 }$$

where the XC part v^xc of the reference mean-field potential has to be subtracted. Equation (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–118].

Rewriting the diagonal elements of the G⁰W⁰ self-energy on the imaginary frequency axis in terms of four-center Coulomb integrals, we obtain

$$\begin{eqnarray} \Sigma_{n\sigma}^{G^0W^0}({\rm i}\omega) &=& \int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \psi_{n\sigma}^\ast({\bf{r}})\Sigma_\sigma^{G^0W^0} ({\bf{r}}, {{\bf{r}}'}, {\rm i}\omega) \psi_{n\sigma}({{\bf{r}}'}) \nonumber \\ &=& -\frac{1}{2\pi} \sum_m \int_{-\infty}^{\infty} {\rm d}\omega^\prime G^0_{m\sigma}({\rm i}\omega+{\rm i}\omega^\prime) (nm,\sigma|W^0({\rm i}\omega^\prime)|mn,\sigma), \end{eqnarray}\noindent \tag{ 29 }$$

where G⁰_mσ(iω) = 1/(iω + _F − _mσ) and

$$\begin{equation} \fl (nm,\sigma|W^0({\rm i}\omega)|mn,\sigma)= \int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} \psi_{n\sigma}^\ast({\bf{r}}) \psi_{m\sigma}({\bf{r}}) W({\bf{r}}, {{\bf{r}}'},{\rm i}\omega) \psi_{m\sigma}^\ast({{\bf{r}}'}) \psi_{n\sigma}({{\bf{r}}'}). \end{equation} \tag{ 30 }$$

This expression can be analytically continued to the real-frequency axis using the Padé approximation or a two-pole model [119].

The G⁰W⁰ approach, as a highly popular choice for quasiparticle excitation calculations, in particular for solids [59–61, 120, 121], is akin to the RPA approach for computing the ground-state correlation energy. Similarly, one can also introduce a second-order perturbation theory for the self-energy [44] consistent with the MP2 correlation energy. Based on the HF non-interacting Green function (equation (24) with HF orbitals), the second-order self-energy can be expressed as

$$\begin{eqnarray} \Sigma^{(2)}(1,2) &=& -{\rm i}\int {\rm d}3 \,{\rm d}4 G^{\rm HF}(1,2)G^{\rm HF}(3,4)G^{\rm HF}(4,3)v(1,3)v(2,4) \nonumber \\ & &+ {\rm i}\int {\rm d}3\,{\rm d}4 G^{\rm HF}(1,4)G^{\rm HF}(4,3)G^{\rm HF}(3,2)v(1,3)v(2,4). \end{eqnarray}\noindent \tag{ 31 }$$

Here for notational simplicity, we have used 1 = (r₁,t₁) and v(1,2) = v(r₁ − r₂)δ(t₁,t₂). After integrating out the internal (spatial and time) coordinates in (31), the final expression for Σ⁽²⁾ (in the frequency domain) within the HF molecular orbital basis is given by

$$\begin{eqnarray} \Sigma^{\rm (2)}_{mn,\sigma} (\omega)& =& \sum_{l}^{\rm occ}\sum_a^{\rm unocc} \sum_{p,\sigma^\prime} (mp,\sigma|la,\sigma^\prime)(pn,\sigma|al,\sigma^\prime) \nonumber \\ & &\times \left[ \frac{\theta(\epsilon_{\rm F}-\epsilon_p)}{\omega+\epsilon_{a\sigma^\prime}-\epsilon_{l\sigma^\prime} - \epsilon_{p\sigma} -{\rm i}\eta} + \frac{\theta(\epsilon_p-\epsilon_{\rm F})} {\omega + \epsilon_{l\sigma^\prime} - \epsilon_{a\sigma^\prime} - \epsilon_{p\sigma} +{\rm i}\eta} \right] \nonumber \\ & & - \sum_{l}^{\rm occ}\sum_a^{\rm unocc} \sum_{p} (mp,\sigma|la,\sigma)(pl,\sigma|an,\sigma) \nonumber \\ & &\times \left[ \frac{\theta(\epsilon_{\rm F}-\epsilon_p)} {\omega+\epsilon_{a\sigma}-\epsilon_{l\sigma} - \epsilon_{p\sigma} -{\rm i}\eta} + \frac{\theta(\epsilon_p-\epsilon_{\rm F})} {\omega + \epsilon_{l\sigma} - \epsilon_{a\sigma} - \epsilon_{p\sigma} + {\rm i}\eta} \right], \end{eqnarray} \tag{ 32 }$$

where θ(x) is the Heaviside step function and η → 0⁺. Again, the two lines in (32) correspond to correlations arising from the second-order direct and second-order exchange interaction, respectively. With this self-energy single-particle excitation energies—here denoted by MP2 quasiparticle energies—can be obtained by adding a correction to the HF orbital energies

$$\begin{equation} \epsilon^{\rm MP2}_{n\sigma} = \epsilon_{n\sigma}^{\rm HF} + \Sigma^{\rm (2)}_{nn,\sigma}(\epsilon^{\rm MP2}_{n\sigma}). \end{equation} \tag{ 33 }$$

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⁰ + G⁰Σ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⁰ and satisfies particle number, momentum and energy conservation [122, 123]. For an in-depth discussion, see [124, 125] and references therein.

In this section we present the basic 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 the other of the above methods; see [76–78, 126–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 (the next sections) in a self-contained way.

The common ingredient of all methods introduced in section 2 (HF, hybrid functionals, MP2, RPA and GW) is the four-orbital Coulomb integrals

$$\begin{equation} (mn,\sigma|ab,\sigma^\prime) = \sum_{ijkl} (ij|kl) c^{i\ast}_{m\sigma} c^{j}_{n\sigma} c^{k\ast}_{a\sigma^\prime} c^{l}_{b\sigma^\prime}, \end{equation} \tag{ 34 }$$

where (ij|kl) are the two-electron integrals in a basis set representation as defined in (14), and cⁱ_mσ are the eigen-coefficients for the molecular orbitals. Computing (and possibly storing) these four-center two-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 have developed the techniques based on multipole expansions of basis functions [132–135]. 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, another 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.

The RI or synonymously the density fitting technique [71–78] amounts to representing pair products of atomic basis functions φ_i(r)φ_j(r) in terms of auxiliary basis functions (ABFs),

$$\begin{equation} \rho_{ij}({\bf{r}}) \equiv \varphi_i({ {\bf{r}}})\varphi_j({ {\bf{r}}}) \approx \tilde{\rho}_{ij}({\bf{r}}) \equiv \sum_{\mu} C_{ij}^{\mu} P_{\mu}({ {\bf{r}}}). \end{equation} \tag{ 35 }$$

μ = 1,2,...,N_aux labels the ABFs {P_μ}, C^μ_ij are the expansion coefficients, and ρ_ij(r) and $\tilde {\rho }_{ij}({\bf{r}})$ here denote pair products of basis functions and their approximate expansion in ABFs. The evaluation of the four-center integrals in (14) then reduces to

$$\begin{equation} (ij|k l) \approx \sum_{\mu\nu} C_{ij}^{\mu} (\mu|\nu) C_{k l}^{\nu}, \end{equation} \tag{ 36 }$$

$$\begin{equation} (\mu|\nu) = V_{\mu\nu}=\int \frac{P_\mu({\bf{r}}) P_\nu({{\bf{r}}'})}{| {\bf{r}}- {{\bf{r}}'}|} {\rm d} {\bf{r}}\, {\rm d} {{\bf{r}}'} . \end{equation} \tag{ 37 }$$

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) four-center integrals reduce to the much cheaper three-center and two-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 (3D) function space, is heavily linearly dependent. Their number scales quadratically with system size, while a non-redundant basis set that expands the same 3D space should scale linearly with system size. For example, the non-interacting response function χ⁰ 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 an 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.

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–78].

The expansion error of basis products in terms of the ABFs (equation (35)) is

$$\begin{equation} \delta \rho_{ij}({\bf{r}})= \tilde{\rho}_{ij}({\bf{r}})-\rho_{ij}({\bf{r}}) = \sum_{\mu} C_{ij}^{\mu} P_{\mu}({\bf{r}}) -\varphi_i({ {\bf{r}}})\varphi_j({ {\bf{r}}}) . \end{equation} \tag{ 38 }$$

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, $\int |\delta \rho _{ij}({\bf{r}})|^2 \mathrm {d} {\bf{r}}$, and yields

$$\begin{equation} C_{ij}^{\mu} = \sum_\nu \langle ij|\nu\rangle S_{\nu\mu}^{-1}, \end{equation} \tag{ 39 }$$

where

$$\begin{equation} \langle ij|\mu\rangle = \int \varphi_i({\bf{r}})\varphi_j({\bf{r}})P_\mu({\bf{r}}){\rm d} {\bf{r}} \end{equation} \tag{ 40 }$$

and

$$\begin{equation} S_{\mu\nu} = \int P_\mu({\bf{r}}) P_\nu({\bf{r}}){\rm d} {\bf{r}} . \end{equation} \tag{ 41 }$$

Combining (36) and (39), one arrives at the following approximation to the four-center Coulomb integrals:

$$\begin{equation} (ij|kl) = \sum_{\mu\nu\mu^\prime\nu^\prime} \langle ij|\mu\rangle S_{\mu\nu}^{-1} V_{\nu\nu^\prime} S_{\nu^\prime\mu^\prime}^{-1} \langle k l |\mu^\prime \rangle . \end{equation} \tag{ 42 }$$

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 four-center integrals themselves,

$$\begin{equation} \delta I_{ij,k l} = (\tilde{\rho}_{ij}|\tilde{\rho}_{k l})-(\rho_{ij}|\rho_{k l}) . \end{equation} \tag{ 43 }$$

As shown by Whitten [72], δI_ij,kl has an upper bound:

$$\begin{eqnarray} &&\fl \delta I_{ij,kl} < (\delta \rho_{ij}|\delta \rho_{ij})^{1/2} (\delta \rho_{k l} | \delta \rho_{k l})^{1/2} + (\tilde{\rho}_{k l}|\delta \rho_{k l})^{1/2} (\delta \rho_{ij}|\delta \rho_{ij})^{1/2} \nonumber\\ &&+ (\tilde{\rho}_{ij}|\delta \rho_{ij})^{1/2} (\delta \rho_{k l} | \delta \rho_{k l})^{1/2} . \end{eqnarray} \tag{ 44 }$$

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–75]

$$\begin{equation} C_{ij}^{\mu} = \sum_{\nu}(ij|\nu)V_{\nu\mu}^{-1} , \end{equation} \tag{ 45 }$$

where

$$\begin{equation} (ij|\nu)=\int\!\!\!\int \frac{\phi_i({\bf{r}})\phi_j({\bf{r}})P_\nu({{\bf{r}}'})}{| {\bf{r}}- {{\bf{r}}'}|} {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} , \end{equation} \tag{ 46 }$$

and the Coulomb matrix V is defined in (37). Combining (36) and (45) one obtains the following decomposition of the four-center ERIs:

$$\begin{equation} (ij|k l) \approx \sum_{\mu\nu} (ij|\mu)V^{-1}_{\mu\nu}(\nu|k l). \end{equation} \tag{ 47 }$$

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–75]. This can be most easily seen by inspecting the error introduced in the self-Coulomb repulsion of the NAO pairs,

$$\begin{eqnarray} &&\delta I_{ij,ij} = (\tilde{\rho}_{ij}|\tilde{\rho}_{ij})- (\rho_{ij}|\rho_{ij}) = 2(\delta\rho_{ij}|\tilde{\rho}_{ij})-(\delta\rho_{ij}|\delta{\rho}_{ij}) . \end{eqnarray} \tag{ 48 }$$

In the RI-V approximation, the first term in the above equation vanishes, and the non-zero contribution comes only from the second order of δρ_ij. This can be readily verified as follows:

$$\begin{eqnarray*} (\delta\rho_{ij}|\tilde{\rho}_{ij}) & = & (\tilde{\rho}_{ij}|\tilde{\rho}_{ij}) - (\rho_{ij}|\tilde{\rho}_{ij}) = \sum_{\mu\nu} C_{ij}^{\mu}V_{\mu\nu}C_{ij}^{\nu} -\sum_{\nu} (ij|\nu)C_{ij}^{\nu} \\ &=& \sum_{\nu} (ij|\nu)C_{ij}^{\nu} - \sum_{\nu} (ij|\nu)C_{ij}^{\nu} = 0, \nonumber \end{eqnarray*}$$

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–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⁴) for HF). 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–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 ABFs 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 HF and hybrid functionals—certainly more than an order of magnitude better than in RI-SVS, creating a competitive alternative for the cases when RI-V is prohibitive. More details go beyond the scope of this paper and will be presented in a forthcoming publication [147].

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):

$$\begin{equation} \Sigma^{\rm x}_{ij,\sigma} = \sum_{k l}\sum_{\mu\nu}(ik|\mu)V_{\mu\nu}^{-1}(\nu|jl)D_{kl} = \sum_{\mu} \sum_{k l} M_{ik}^{\mu}M_{jl}^{\mu}D_{k l}, \end{equation} \tag{ 49 }$$

where

$$\begin{equation} M_{ik}^{\mu} = \sum_{\nu} (ik|\nu)V_{\nu\mu}^{-1/2} =\sum_{\nu} C_{ik}^{\nu} V_{\nu\mu}^{1/2} . \end{equation} \tag{ 50 }$$

Σ^x_ij,σ must thus be recomputed for each iteration within the self-consistent field (scf) loop. The required floating point operations scale as N³_b·N_aux, where N_b is the number of regular single-particle basis functions. The transformation matrix M^μ_ik, on the other hand, is constructed only once (prior to the scf loop), requiring N²_b·N²_aux operations.

The numerical efficiency can be further improved by inserting the expression for the density matrix (13) into (49):

$$\begin{equation} \Sigma^{\rm x}_{ij,\sigma}=\sum_{n}^{\rm occ}\sum_{\mu}\left(\sum_{k}M_{ik}^{\mu} c_{n \sigma}^k\right)\left(\sum_{l} M_{jl}^{\mu}c_{n \sigma}^l\right) = \sum_{n}^{\rm occ}\sum_{\mu}B_{in\sigma}^{\mu}B_{jn\sigma }^{\mu} . \end{equation} \tag{ 51 }$$

The formal scaling in (51) is now N_occ·N_b²·N_aux, with N_occ being the number of occupied orbitals, and is thus improved by a factor N_b/N_occ (typically 5–10). Once the exact-exchange matrix is obtained, the EX follows through

$$\begin{equation} E_{\rm x}^{\rm HF} = -\frac{1}{2} \sum_{ij,\sigma} \Sigma^{\rm x}_{ij,\sigma} D_{ij,\sigma} . \end{equation} \tag{ 52 }$$

For a variety of physical problems, combining 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:

$$\begin{equation} E^{\rm hyb}_{\rm xc}=E^{\rm GGA}_{\rm xc} + \alpha (E^{\rm HF}_x - E^{\rm GGA}_{\rm x}) . \end{equation} \tag{ 53 }$$

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 the HF exchange is implemented, it is straightforward to also perform hybrid functional calculations.

To compute the MP2 correlation energy in (16) and the MP2 self-energy in (32) using the RI technique, the MO-based four-orbital two-electron Coulomb integrals are decomposed as follows:

$$\begin{equation} (ma,\sigma|nb,\sigma^\prime)=\sum_{\mu}O^{\mu}_{ma,\sigma}O^{\mu}_{nb,\sigma^\prime} . \end{equation} \tag{ 54 }$$

The three-orbital integrals can be evaluated by

$$\begin{equation} O^{\mu}_{ma,\sigma}=\sum_{ij}M^{\mu}_{ij}c_{m\sigma}^{i\ast}c_{a\sigma}^j \end{equation} \tag{ 55 }$$

following (34), (47) and (50). Plugging (54) into (16), one obtains the RI-V version of the MP2 correlation energy

$$\begin{eqnarray} E^{(2)}_0 & = & \frac{1}{2} \sum_{mn}^{\rm occ}\sum_{ab}^{\rm unocc}\sum_{\sigma,\sigma^\prime} \left(\sum_{\mu}O^{\mu}_{ma,\sigma}O^{\mu}_{nb,\sigma^\prime}\right) \nonumber \\ & & \times \left[ \frac{ \left(\sum_{\mu}O^{\mu}_{am,\sigma}O^{\mu}_{bn,\sigma^\prime}\right)- \left(\sum_{\mu}O^{\mu}_{bm,\sigma}O^{\mu}_{an,\sigma}\right)\delta_{\sigma\sigma^\prime}} {\epsilon_{m\sigma} + \epsilon_{n\sigma^\prime} -\epsilon_{a\sigma}-\epsilon_{b\sigma^\prime}} \right] . \end{eqnarray} \tag{ 56 }$$

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²_b·N_aux and the summation in (56) as N²_occ·(N_b − N_occ)²·N_aux. Like in HF, the scaling exponent ${\mathcal{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 analogous to that of the correlation energy. In our implementation, we first calculate the MP2 self-energy on the imaginary frequency axis, Σ⁽²⁾_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 be used to cross-check each other to guarantee the reliability of the final results.

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 $\left [\chi ^0(\mathrm {i}\omega)v\right ]$. Combining (19) and (54) we obtain

$$\begin{equation} {\rm Tr}\left[\chi^0({\rm i}\omega)v\right]= \sum_{\mu} \sum_\sigma \sum_m^{\rm occ}\sum_a^{\rm unocc} \frac{O^\mu_{ma,\sigma}O^\mu_{am,\sigma}} {{\rm i}\omega-\epsilon_{a\sigma}+\epsilon_{m\sigma}} + {\rm c.c.} , \end{equation} \tag{ 57 }$$

where O^μ_ma,σ is given by (55). Next we introduce an auxiliary quantity Π(iω):

$$\begin{equation} \Pi({\rm i}\omega)_{\mu\nu} = \sum_\sigma \sum_m^{\rm occ}\sum_a^{\rm unocc} \frac{O^\mu_{ma,\sigma}O^\nu_{am,\sigma}} {{\rm i}\omega-\epsilon_{a\sigma}+\epsilon_{m\sigma}} + {\rm c.c.} , \end{equation} \tag{ 58 }$$

which allows us to write

$$\begin{equation} {\rm Tr}\left[\chi^0({\rm i}\omega)v\right]={\rm Tr}\left[v^{1/2}\chi^0({\rm i}\omega)v^{1/2}\right] = {\rm Tr}\left[\Pi({\rm i}\omega)\right] . \end{equation} \tag{ 59 }$$

Thus the matrix Π can be regarded as the matrix representation of the composite quantity v^1/2χ⁰(iω)v^1/2 using the ABFs. It is then easy to see that the RPA correlation energy (18) can be computed as

$$\begin{eqnarray} E_{\rm c}^{\rm RPA} &= & \frac{1}{2\pi}\int_0^\infty {\rm d}\omega {\rm Tr} \left[{\rm ln}\left(1-\Pi({\rm i}\omega)\right) + \Pi({\rm i}\omega)\right] \nonumber \\ &=& \frac{1}{2\pi} \int_0^\infty {\rm d}\omega \left\{ \ln\left[{\rm det}(1-\Pi({\rm i}\omega))\right] + {\rm Tr}\left[ \Pi({\rm i}\omega) \right] \right\} \end{eqnarray} \tag{ 60 }$$

using the general property Tr $\left [ \ln (A)\right ]=\ln \left [{\mathrm { det}}(A)\right ]$ 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 C for further details) with typically 20–40 frequency points. For a fixed frequency grid size, the number of required operations is proportional to N_occ·N_unocc·N²_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 ${\mathcal{O}}(N^4)$ scaling. An ${\mathcal{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 Π(iω) = v^1/2χ₀(iω)v^1/2, we obtain

$$\begin{eqnarray*} W_{\mu\nu}({\rm i}\omega) &=& \int\!\!\!\int {\rm d} {\bf{r}} \,{\rm d} {{\bf{r}}'} P_\mu^\ast({\bf{r}})W({\bf{r}}, {{\bf{r}}'},{\rm i}\omega)P_\nu({\bf{r}}') \\ &=& \sum_{\mu^\prime \nu^\prime} V^{1/2}_{\mu\mu^\prime} [1-\Pi({\rm i}\omega)]^{-1}_{\mu^\prime \nu^\prime} V^{1/2}_{\nu^\prime\nu} . \end{eqnarray*}$$

To apply the RI-decomposition to (30), we expand

$$\begin{equation} \psi_{n\sigma}^\ast({\bf{r}}) \psi_{m\sigma}({\bf{r}})=\sum_{ij}\sum_{\mu} P_{\mu}({\bf{r}}) C_{ij}^{\mu}c_{n\sigma}^{i\ast} c_{m\sigma}^j , \end{equation} \tag{ 61 }$$

where (10) and (35) are used. Combining (30), (50), (55) and (61) gives

$$\begin{equation} (nm,\sigma|W^0({\rm i}\omega)|mn,\sigma)= \sum_\sigma \sum_{\mu\nu} O_{nm,\sigma}^\mu \left[1-\Pi({\rm i}\omega)\right]^{-1}_{\mu\nu} O_{mn,\sigma}^\nu . \end{equation} \tag{ 62 }$$

Inserting (62) into (29), one finally arrives at the RI version of the G⁰W⁰ self-energy

$$\begin{eqnarray} \Sigma_{n\sigma}^{G^0W^0}({\rm i}\omega) &= & -\frac{1}{2\pi} \sum_m \int_{-\infty}^{\infty} {\rm d}\omega^\prime \frac{1}{{\rm i}\omega+{\rm i}\omega^\prime+\epsilon_{\rm F}-\epsilon_{m\sigma}} \nonumber \\ & & \times \sum_{\mu\nu} O_{nm,\sigma}^\mu \left[1-\Pi({\rm i}\omega)\right]^{-1}_{\mu\nu} O_{mn,\sigma}^\nu . \end{eqnarray} \tag{ 63 }$$

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).

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 in 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

$$\begin{equation} \varphi_i({\bf{r}}) = \frac{u_{s(a) l \kappa}(r)}{r}Y_{lm}(\hat{ {\bf{r}}}_a) . \end{equation} \tag{ 64 }$$

u_s(a)lκ is a radial function centered at atom a, and $Y_{lm}(\hat { {\bf{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 $\hat { {\bf{r}}}_a=({\bf{r}}- {\bf{R}}_a)/| {\bf{r}}- {\bf{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 paper.

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 per 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 also use, e.g., LDA minimal basis functions for calculations using other functionals if their minimal basis functions are not available. Below we discuss this possibility for HF.

To give one specific example, consider the nitrogen atom (this case and more are spelled out in detail in table 1 of [64]). There are five minimal basis functions, of 1s, 2s and 2p orbital character, respectively. In shorthand notation, we denote the number of radial functions for given angular momenta by (2s1p) (two s-type radial functions and 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 (4s3p2d1f1g) in tier 2, 55 (5s4p3d2f1g) in tier 3 and 80 (6s5p4d3f2g) in tier 4.

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, 127, 152]. 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 similar 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

$$\begin{equation} P_{\mu}({\bf{r}}) = \frac{\xi_{s(a) l \kappa }(r)}{r}Y_{lm}(\hat{ {\bf{r}}}_a) \end{equation} \tag{ 65 }$$

just like for the one-particle NAO basis functions in (64), but of course with different radial functions. To distinguish the ABFs from the NAO basis functions, we denote the radial functions of the ABFs by ξ_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 ABFs follows from the definition of the orbital basis set. At each level of NAO basis, one can generate a corresponding ABF set, denoted by 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 $\{\tilde {\xi }_{s l\kappa }(r)=u_{s k_1l_1}(r)u_{s k_2l_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_sk1l1 and u_sk2l2, i.e. |l₁ − l₂| ⩽ l ⩽ |l₁ + l₂|.
• (ii)  
  Even for relatively small orbital basis sets, the number of resulting auxiliary radial functions $\{\tilde {\xi }_{sl \kappa }(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 by the spherical harmonics $Y_{lm}(\hat { {\bf{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 dependences 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 paper (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⁻² or 10⁻³ 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⁻⁴ or 10⁻⁵, 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.

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 very important, but since many pieces of our eventual implementation exist in the literature, here we give only a brief summary; see the 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 (appendix A.1). The required three-center integrals

$$\begin{equation} (ij|\mu) = \int \phi_i({\bf{r}})\phi_j({\bf{r}})Q_\mu({\bf{r}}){\rm d} {\bf{r}} \end{equation} \tag{ 66 }$$

are carried out by standard overlapping atom-centered grids as used in many quantum-chemical codes for the XC matrix in DFT [64, 65, 154, 155]; see appendix A.2. The same strategy works for two-center integrals

$$\begin{equation} (\mu|\nu) = \int P_\mu({\bf{r}})Q_\nu({\bf{r}}){\rm d} {\bf{r}}. \end{equation} \tag{ 67 }$$

As an alternative, we have also implemented two-center integrals following the ideas developed by Talman [133, 134], which are described in appendices A.3 and A.4. In summary, we thus have accurate matrix elements at hand that are used for the remainder of this work.

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₂ 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⁻², and the global SVD threshold ε^svd is set to 10⁻⁴, 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: firstly, the quality of the ABF expansion improves as the number of ABFs increases. Secondly, 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 behavior 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 3D integrations, the integration weight r² dr is small).

Zoom In Zoom Out Reset image size

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 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 in particular 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 larger effect on RI-SVS than RI-V, underscoring the desired variational properties of RI-V [71–75].

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 the results presented in this section correspond to the cc-pVQZ basis set of Dunning et al [149, 150, 157]. The convergence behavior with respect to NAO/GTO single-particle basis set size will be the topic of the next section.

We first check the quality of the ABF prescription for light (first and second row) elements. We again choose N₂ as a first illustrating example. Table 2 presents RI-HF total and binding energy errors for the N₂ 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 of thresholding parameters ε^orth and ε^svd (see section 4.2). All binding energy errors are obtained after a counterpoise correction [81] to remove any possible basis set superposition errors (BSSE) (see also the lowest panel of figure 2).

Zoom In Zoom Out Reset image size

For N₂ 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₂. 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 total-energy 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⁻², ε^svd = 10⁻⁴) (standard accuracy) and (ε^orth = 10⁻³, ε^svd = 10⁻⁵) (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₂ (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 figures 3 and 4 the non-relativistic RI-V HF and MP2 total energy errors and atomization energy errors for 20 molecules are shown. In all cases the total energy error is below 1 meV per atom, demonstrating that the meV accuracy 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 versus 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^{{\rm ABF}\hbox{-}{\rm max}}$) to be just 1 higher than that of the single-particle atomic orbitals ($l^{{\rm AO}\hbox{-}{\rm max}}$).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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₂ as a function of the bond length. Again the cc-pVQZ basis and the NWChem reference values are used here.

Zoom In Zoom Out Reset image size

Using the thresholding parameters (ε^orth = 10⁻², ε^svd = 10⁻⁴), the RI-V HF total energy error can be as large as ∼15 meV per atom for copper, in contrast to the <1 meV per 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 per atom total-energy accuracy can be achieved at ε^orth = 10⁻⁴, 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 per atom along a large range of bonding distances, regardless of the choice of thresholding parameters. Also, increasing the highest angular momentum for ABFs beyond $l^{{\rm AO}\hbox{-}{\rm 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₂ 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⁻⁵, ε^svd = 10⁻⁵). 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.

Zoom In Zoom Out Reset image size

From figure 6 one can see that with thresholding parameters that are sufficient for light elements (ε^orth = 10⁻², ε^svd = 10⁻⁴), the error in the total energy is even bigger—one order of magnitude larger than in the case of Cu₂. However, by going to tighter and tighter on-site ABF thresholding parameter ε^orth, the total-energy error can again be reduced to the meV per atom level. Similar to the Cu₂ 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 reliable results for the RI-V approximation. For light elements, one can readily get meV per atom accuracy in total energies and sub-meV per atom accuracy in binding energies, for a wide range of thresholding parameters. For heavy elements, meV per atom accuracy requires tighter thresholding parameters, particularly 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⁻² in the following. For heavier elements (Z > 18), we resort to ε^orth = 10⁻⁴, which yields negligibly small errors in total energies even for heavy elements, and ε^orth = 10⁻³ for elements in between. The additional, system-wide SVD threshold ε^svd is set to 10⁻⁴ or tighter for the remainder of this paper. As shown above, its accuracy implications are then negligible as well.

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 functional-independent optimized basis sets (tiers). Thus, besides the discussion of the convergence behavior of the generic optimized tier basis sets, here we will also address the influence of the choice of the minimal basis which is in practice generated by certain atomic solvers. 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 behavior 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 were generated by the HF atomic solver. Unfortunately, at the moment a 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 tier 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 [102, 160], 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 above-mentioned problem is alleviated to some extent.

5.1.1. Light elements

We first check the convergence behavior 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₂ 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 behavior of the LDA total energy of N₂ 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 CBS limit [162]³.

Zoom In Zoom Out Reset image size

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. If the true HF minimal basis were used, we should expect an even faster basis convergence of the HF total energy, similar to the LDA total-energy convergence behavior 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₂ 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 tiers 2 and 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 from 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-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).

Zoom In Zoom Out Reset image size

Figure 9 demonstrates the same behavior for a different test case, the binding energy of the water dimer using HF (left) and the PBE0 [7, 30] 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₂O fragments with fixed geometry as in the dimer, not to fully relaxed H₂O monomers. This is sufficient for the purpose of the basis convergence test here. Detailed geometrical information for the 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₂, 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 the exact exchange is included.

Zoom In Zoom Out Reset image size

In practice, HF calculations at the tier 2 level of our NAO basis sets yield accurate results for light elements. Counterpoise corrections help us 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).

5.1.2. 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₂, Ag₂ and Au₂ at a fixed binding distance, d = 2.5 Å. (The experimental binding distances are 2.22 Å [163], 2.53 Å [164, 165] and 2.47 Å [163], respectively.) Again, we find that KLI-derived minimal basis sets are noticeably better converged (lower total energies) than LDA-derived minimal basis sets. In contrast to N₂, 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).

Zoom In Zoom Out Reset image size

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 co-workers [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 errors 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 chemistry 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 [166, 167]. The agreement between our all-electron approach and the GTO-ECP one is pretty decent for Ag₂ and Au₂, whereas a larger discrepancy of ∼0.03 eV is observed for Cu₂. For this particular case, we suspect that 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.

In the implementation described here, the 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 [38, 86, 90, 168–174]. Much experience has been gained in the quantum chemistry community to construct Gaussian basis sets for correlated calculations [157, 175], but for NAOs this is not the 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 [38, 86, 90, 168–171], 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₂. 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 topic of 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 (appendix C), which ensures a high accuracy for the systems studied here.

Zoom In Zoom Out Reset image size

5.2.1. Binding energies

As illustrating examples for light elements, in figure 12 the BSSE-corrected MP2 and RPA binding energies for N₂ 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 those 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₂ [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₂ 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.

Zoom In Zoom Out Reset image size

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 as '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 and −206.5 meV for the MP2 and RPA@HF binding energies, comparable to the quality of the cc-pV6Z basis that yields −221.1 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 behavior 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 reduced or increased the onset radii of this confining potential? The answer to this question is illustrated in figure 13 where the basis convergence behaviors for N₂ 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 behavior for N₂ or (H₂O)₂. Upon reducing it, noticeable changes in 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 behavior described above. This finding holds in general for covalent and hydrogen bonds. In practice, the onset radius 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 Å) in our experience. Further details of this can be found in [176].

Zoom In Zoom Out Reset image size

We next illustrate the NAO basis convergence for heavy elements, using Au₂ as an example. In figure 14 the MP2 binding energy for the Au₂ dimer as a function of the bond length is plotted for different NAO basis sets. Relativity is again treated at the scaled ZORA level [64, 158]. The binding curves shown here demonstrate that the same qualitative convergence behavior 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 tiers 2 and 4 yield essentially converged results. For comparison, we show a completely independent (NWChem calculations) curve with Gaussian 'aug-cc-pV5Z-PP' basis sets [166, 167]. The resulting binding energy curve yields rather close agreement with our all-electron, NAO basis set results.

Zoom In Zoom Out Reset image size

5.2.2. 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⁰W⁰@HF and MP2 quasiparticle HOMO levels for N₂ 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 an MP2 total-energy 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 still considerable, ∼0.2–0.3 eV for G⁰W⁰ and ∼0.1–0.2 eV for MP2-QP. These errors are brought down to ∼0.1 eV for G⁰W⁰ and ∼0.06 eV for MP2-QP by going to a pure tier 4 NAO basis set. 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 [38, 86, 90, 168–171]. Therefore 'tier 4 + a5Z-d' basis sets were used in the benchmark calculations presented in section 6.

Zoom In Zoom Out Reset image size

5.2.3. G⁰W⁰ calculations for the benzene molecule

The G⁰W⁰ method has been used to calculate the single-particle properties for various small- and medium-size molecules with considerable success [128, 153, 186, 188, 191–200]. The dependence of G⁰W⁰ on the starting point has also been looked at in the past [60, 61, 195, 200]. Here, using a selected collection of atoms and molecules from the G2 ion test set for IEs [79], we aim to establish the overall performance of G⁰W⁰ 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 vertical IEs. With vertical IEs we denote IEs at fixed geomerty, i.e. no structural relaxation after excitation. This quantity is directly comparable to the GW and MP2-QP 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 performed 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⁰W⁰@HF and G⁰W⁰@PBE0 for this database. The actual IE values are presented in appendix B. Compared to HF, one can see that the deviations from the experimental values are much smaller in G⁰W⁰ and MP2-QP. On average HF tends to overestimate IEs. This trend is corrected by MP2 and G⁰W⁰. Using the same HF reference, the magnitude of the correction is smaller in G⁰W⁰ 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⁰W⁰ method, G⁰W⁰ based on HF gives too large IEs on average, whereas G⁰W⁰ based on PBE does the opposite. Consequently, G⁰W⁰ 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 yield very similar results, implying that for the light elements reported here orbital relaxation effects are not significant⁴. Table 4 summarizes the error statistics. For the subset of atoms and molecules, MP2-QP gives the smallest mean error (ME) and G⁰W⁰@PBE0 the smallest mean absolute error (MAE).

Zoom In Zoom Out Reset image size

The 55 atomization energies from the original G2-I set (table 3 of [201]) 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 [201], i.e. all-electron MP2 optimization with a Gaussian 6-31G* basis set. The reference data, taken from [202], are 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 [203] and Paier et al [90], respectively, both obtained with GTO basis sets and extrapolated to the CBS 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 behavior of the standard RPA scheme have been proposed. These comprise contributions from single excitations [40] and second-order screened exchange [90, 110].

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 a 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 [204], which we therefore adopt as the reference here. RPA-type calculations for S22 have recently been reported [40, 98, 205]. 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 are presented. 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 an MAE of 39 meV and an 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 behavior in [40] and proposed a simple solution to overcome this problem.

To reduce the (formally) 6D Coulomb integrations to 3D ones, we first solve Poisson's equation for each P_μ(r') (classical electrostatics). We define Q_μ(r) as

$$\begin{equation} Q_{\mu}({\bf{r}}) = \int {\rm d}{ {{\bf{r}}'}} v({\bf{r}}- {{\bf{r}}'})P_\mu({{\bf{r}}'}). \end{equation} \tag{ A.1 }$$

Using the Laplace expansion of the Coulomb potential v(r − r') = 1/|r-r'|,

$$\begin{equation} v({\bf{r}}- {{\bf{r}}'})=\sum_{lm}\frac{4\pi}{2l+1}Y_{lm}(\hat{\bf{r}})Y^\ast_{lm}(\hat{\bf{r}}') \times \frac{r_{<}^l}{r_{>}^{l+1}}, \end{equation} \tag{ A.2 }$$

where r_< = min(r,r') and r_> = max(r,r'), the integration of the angular part in (A.1) can be done analytically, yielding

$$\begin{equation} Q_{\mu}({\bf{r}}) = \frac{\alpha_{al\kappa}(r)}{r}Y_{lm}(\hat{ {\bf{r}}}_a). \end{equation} \tag{ A.3 }$$

The radial part α_alκ(r) is given by a simple 1D (numerical) integration [64, 65, 209]

$$\begin{equation} \alpha_{al\kappa}(r) = \frac{4\pi}{2l+1} \left[ \left(\int _0^{r} {\rm d} r_< \xi_{al\kappa}(r_<) \frac{r_<^{l+1}}{r^{l}}\right) + \left(\int _r^{\infty} {\rm d} r_> \xi_{al\kappa}(r_>) \frac{r^{l+1}}{r_>^{l}}\right) \right]. \end{equation} \tag{ A.4 }$$

Thus, the three-center and two-center Coulomb integrals (37) and (45) reduce to the 3D integrals in (66) and (67).

The three-center and two-center integrals in this work are computed by grid-based integrations using overlapping, atom-centered spherical grids and the same technology that is used in many quantum-chemical applications for the XC matrix in the DFT [65, 154]. The integration grid points r = r(a,s,t) are uniquely specified by the atomic center a, the radial shell number s and the 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₁, a₂, a₃, the three-center integrals in (66) can then be discretized as

$$\begin{equation} (ij|\mu) = \sum_{a={a_1,a_2,a_3}} \sum_{s,t}p_3(a, {\bf{r}}) w(s,t) \phi_i({\bf{r}})\phi_j({\bf{r}})Q_\mu({\bf{r}}), \end{equation} \tag{ A.5 }$$

where p₃(a,r) is a three-center partition function that satisfies

$$\begin{eqnarray*} &&\sum_{a=a_1,a_2,a_3}p_3(a, {\bf{r}})=1 \end{eqnarray*}$$

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.

As mentioned in section 4.3, two-center integrals of numeric atom-centered basis functions like V_μν in (37) and S_μν in (41) can be efficiently calculated in Fourier space as described by Talman [132, 210]. This and the following (appendix A.4) subsections give the details of our implementation. We first describe the general procedure here, and the central ingredient of our implementation, the logarithmic spherical Bessel transform (logarithmic SBT (logSBT)), will be presented in the next subsection.

As is well known, the overlap of two functions f(r) and g(r) can be expressed in Fourier space as

$$\begin{equation} \int f({\bf{r}}) g({\bf{r}} - \boldsymbol{R}) \,{\rm d}{\bf{r}} = \int {\skew3\tilde{f}}(\boldsymbol{k}) {\skew3\tilde{g}}(-\boldsymbol{k})\, {\rm e}^{{\rm i}\boldsymbol{k}\cdot\boldsymbol{R}} \,{\rm d}\boldsymbol{k}. \end{equation} \tag{ A.6 }$$

The Fourier transform ${\skew3\tilde{f}}(\boldsymbol{k})$ in (A.6) of an atomic function $f({\bf{r}}) = f(r) Y_{lm}(\hat {\bf{r}})$ has the same angular momentum in real and Fourier space for symmetry reasons

$$\begin{equation} {\skew3\tilde{f}}(\boldsymbol{k}) = (2\pi)^{-{3}/{2}} \int f({\bf{r}}) {\rm e}^{-{\rm i}\boldsymbol{k}\cdot{\bf{r}}} \,\,{\rm d} {\bf{r}} = {\rm i}^{-l} {\skew3\tilde{f}}(k) Y_{lm}(\hat{\boldsymbol{k}}). \end{equation} \tag{ A.7 }$$

The radial part ${\skew3\tilde{f}}(k)$ is given by the SBT of f(r)

$$\begin{equation} {\skew3\tilde{f}}(k) = \sqrt{\frac{2}{\pi}} \int_0^\infty j_{l}(kr) f(r)\,r^2 \,{\rm d} r. \end{equation} \tag{ A.8 }$$

If f(r) is tabulated on a logarithmic grid, its SBT can be calculated efficiently on an equivalent logarithmic grid using the fast logSBT algorithm [211–213] as described in the next section.

The 3D integral (A.6) can be separated by expanding the plane wave e^ik·R in spherical harmonics and spherical Bessel functions

$$\begin{equation} {\rm e}^{{\rm i} \boldsymbol{k}\cdot\boldsymbol{R}} = (2\pi)^{{3}/{2}} \sum_{L=0}^\infty {\rm i}^L \sqrt{\frac{2}{\pi}} j_L(kR) \sum_M Y_{LM}(\hat{\boldsymbol{k}}) Y_{LM}^*(\hat{\boldsymbol{R}}). \end{equation} \tag{ A.9 }$$

The separation yields

$$\begin{equation} \int f({\bf{r}}) g({\bf{r}} - \boldsymbol{R})\, \,{\rm d}{\bf{r}} = (2\pi)^{{3}/{2}} \sum_L {\rm i}^{-l+l'+L} I_L(R) A_L(\hat{\boldsymbol{R}}) \end{equation} \tag{ A.10 }$$

with a radial integral

$$\begin{equation} I_L(R) = \sqrt{\frac{2}{\pi}} \int_0^\infty j_L(kR) {\skew3\tilde{f}}(k) {\skew3\tilde{g}}(k)\, k^2 \,{\rm d} k \end{equation} \tag{ A.11 }$$

and an angular integral

$$\begin{equation} A_L(\hat{\boldsymbol{R}}) = \sum_M Y_{LM}^*(\hat{\boldsymbol{R}})\, C\bigl(lm; l'm'; LM\bigr). \end{equation} \tag{ A.12 }$$

The triple-Y integrals

$$\begin{equation} C\bigl(lm; l'm'; LM\bigr) := \int_{\Omega} Y_{lm}(\hat{\boldsymbol{k}}) Y_{l'm'}(\hat{\boldsymbol{k}}) Y_{LM}(\hat{\boldsymbol{k}}) \, \,{\rm d} \hat{\boldsymbol{k}} \end{equation} \tag{ A.13 }$$

in (A.12) can be calculated efficiently using recursion formulae [210]. They are non-zero 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 a given distance R, the radial integrals I_L(R) in (A.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 (A.11) as an SBT of $\tilde P(k) = {\skew3\tilde{f}}(k) {\skew3\tilde{g}}(k)$. By applying logSBT, and interpolating for all needed distances R, one can obtain all the integrals at tight-binding cost, meaning that an efficient recursion formula (together with spline evaluations) can be used instead of evaluating each integral numerically.

Coulomb interactions of atomic functions can be calculated with comparable ease where the integrand in (A.11) is multiplied with the Coulomb kernel 4π/k²

$$\begin{equation} V_L(R) = 4\pi \sqrt{\frac{2}{\pi}} \int_0^\infty j_L(kR) \frac{{\skew3\tilde{f}}(k) {\skew3\tilde{g}}(k)}{k^2}\,k^2\,{\rm d} k. \end{equation} \tag{ A.14 }$$

The Coulomb interaction generally does not vanish even if the two charge densities do not overlap. However, it has a simple multipolar behavior and explicit integration of (A.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 SBT $\tilde P(k) = {\skew3\tilde{f}}(k){\skew3\tilde{g}}(k)$. Therefore, it only depends on the multipole moment of P, which is encoded in its limiting behavior 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

$$\begin{equation} V_{l+l'}(R) = \frac{4\pi}{R^{l+l'+1}} \sqrt{\frac{2}{\pi}} \frac{(2l+2l'-1)!!}{(2l-1)!!\, (2l'-1)!!}\, p_f p_g \end{equation} \tag{ A.15 }$$

with (2n − 1)!! = 1·3···(2n − 1) and with multipole moments

$$\begin{equation} p_f = \frac{1}{2l+1} \int_0^\infty r^{2+l} f(r)\,\,{\rm d} r. \end{equation} \tag{ A.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 et al in [214], 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.

This section describes our implementation of the logSBT algorithm. For a more extensive description, see Talman [213] and Hamilton [212].

The SBT as defined in (A.8) can be written as the integral over a kernel $\mathcal J(kr)$ and a right-hand-side $\mathcal F(r)$:

$$\begin{equation} {\skew3\tilde{f}}(k) = k^{-\alpha} \int_0^\infty \frac{\,{\rm d} r}{r} \underbrace{\sqrt{\frac{2}{\pi}} j_l(kr) (kr)^\alpha}_{=: \mathcal J(kr)} \, \underbrace{r^{3-\alpha} f(r)}_{=: \mathcal F(r)}. \end{equation} \tag{ A.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 $\mathcal J(\kappa + \rho)$ and (A.17) turns into a convolution, which can be efficiently calculated using fast Fourier transforms (FFTs). Please note that dr/r = dρ in (A.17).

As pointed out by Hamilton [212], this procedure is exact if both $\mathcal F(\rho) = \mathrm {e}^{(3-\alpha)\rho } f(\mathrm {e}^{\rho })$ and the corresponding SBT term $\tilde {\mathcal F}(\kappa) = \mathrm {e}^{\alpha \kappa } {\skew3\tilde{f}}(\mathrm {e}^{\kappa })$ 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 $\mathcal F(\rho)$ and $\tilde {\mathcal F}(\kappa)$ smoothly drop to zero on both ends, which therefore can safely be connected.

Unfortunately, the scaling factor k^−α turns out to be quite problematic. By the design of the algorithm, the absolute error of $\tilde F(\kappa) = k^{\alpha } f(k)$ before the final scaling is typically of the order of machine precision, i. e. about 10⁻¹⁵. 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 [213] 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 $\mathcal F(\rho)$ 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^\alpha {\skew3\tilde{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₀(e^τ) into a smooth part proportional to erfc(τ/Δτ₀) and a properly decaying rest. We use the first part for a trapezoidal rule and the second part for the log-periodic algorithm. Fortunately, 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 [213], 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 behavior.