Separable Resolution-of-the-Identity with All-Electron Gaussian Bases: Application to Cubic-scaling RPA

TL;DR

This paper introduces a cubic-scaling separable resolution-of-the-identity method for all-electron Gaussian basis calculations, enabling efficient RPA correlation energy computations with maintained accuracy.

Contribution

It develops a novel density fitting scheme using quadratures over real-space points that preserves atomic orbital use and scales cubically for all-electron RPA calculations.

Findings

Accurate Fock exchange and MP2 energies demonstrated.

Cubic scaling in operations for RPA correlation energies.

Small crossover point with standard RI-RPA implementation.

Abstract

We explore a separable resolution-of-the-identity formalism built on quadratures over limited sets of real-space points designed for all-electron calculations. Our implementation preserves in particular the use of common atomic orbitals and their related auxiliary basis sets. The set up of the present density fitting scheme, i.e. the calculation of the system specific quadrature weights, scales cubically with respect to the system size. Extensive accuracy tests are presented for the Fock exchange and MP2 correlation energies. We finally demonstrate random phase approximation (RPA) correlation energy calculations with a scaling that is cubic in terms of operations, quadratic in memory, with a small crossover with respect to our standard RI-RPA implementation.

Figures

Tables

Table 1. Table 1: RPA correlation energy (Hartree) for the benzene molecule. The various RI schemes are compared. Calculations are performed at the (spherical) cc-pVTZ level with the corresponding (cartesian) cc-pVTZ-RI auxiliary basis. The number of imaginary frequency in the RPA integration is given by n ω 𝜔 \omega . The number of times in the Laplace transform (LT) approach is given by n τ 𝜏 \tau .

n$\omega$	RI-SVS	RI-V	RI-RS	RI-RS+LT	RI-RS+LT
				(n$\tau$=1.5$\times$n$\omega$)	(n$\tau$=2$\times$n$\omega$)
6	-1.25169968	-1.25148948	-1.25143651	-1.25148267	-1.25143097
8	-1.25171657	-1.25150621	-1.25145324	-1.25145455	-1.25145327
10	-1.25171791	-1.25150754	-1.25145456	-1.25145456	-1.25145456
12	-1.25171791	-1.25150754	-1.25145456	-1.25145456	-1.25145456

Equations

$$\begin{array}{rl}\varphi (\mathbf{r}){\varphi }^{\prime }(\mathbf{r})\simeq & \sum _{\beta }{\mathcal{F}}_{\beta }(\varphi {\varphi }^{\prime })\beta (\mathbf{r})\\ \doteqdot & \mathcal{F}(\varphi {\varphi }^{\prime };\mathbf{r})\end{array}$$

$${\mathcal{F}}_{\beta }^V(\varphi {\varphi }^{\prime })=\sum _{{\beta }^{\prime }}[V^{-1}{]}_{\beta {\beta }^{\prime }}({\beta }^{\prime }|\varphi {\varphi }^{\prime })$$

$${\mathcal{F}}_{\beta }^{SVS}(\varphi {\varphi }^{\prime })=\sum _{{\beta }^{\prime }}[S^{-1}{]}_{\beta {\beta }^{\prime }}⟨{\beta }^{\prime }|\varphi {\varphi }^{\prime }⟩$$

$$(\beta |\varphi {\varphi }^{\prime })$$

$$⟨\beta |\varphi {\varphi }^{\prime }⟩$$

$$\begin{array}{rl}{\mathcal{F}}_{\beta }(\varphi {\varphi }^{\prime })& =\sum _k[M{]}_{\beta k}⟨f_k|\varphi {\varphi }^{\prime }⟩\\ & =\sum _k[M{]}_{\beta k}⟨f_k|\varphi ⟩⟨f_k|{\varphi }^{\prime }⟩\end{array}$$

$$\begin{array}{rl}{\mathcal{F}}_{\beta }^{RS}(\varphi {\varphi }^{\prime })& =\sum _k[M{]}_{\beta k}⟨\delta (\mathbf{r}-{\mathbf{r}}_k)|\varphi {\varphi }^{\prime }⟩\\ & =\sum _k[M{]}_{\beta k}\varphi ({\mathbf{r}}_k){\varphi }^{\prime }({\mathbf{r}}_k)\end{array}$$

$$\operatorname*{argmin}_{M}\sum_{\rho,\beta}\Big{(}\mathcal{F}_{\beta}^{RS}(\rho)-\mathcal{F}_{\beta}^{V}(\rho)\Big{)}^{2}$$

$$\operatorname*{argmin}_{M}\Big{|}\Big{|}M\cdot D-F\Big{|}\Big{|}$$

$$M=F\cdot D^{†}\cdot (D\cdot D^{†}{)}^{-1}$$

$$(D\cdot D^{†}{)}^{-1}\simeq d\cdot (\tilde{D}\cdot {\tilde{D}}^{†}+ϵ\text{mathds}I{)}^{-1}\cdot d$$

$$M=F\cdot {\tilde{D}}^{†}\cdot (\tilde{D}\cdot {\tilde{D}}^{†}+ϵ\text{mathds}I{)}^{-1}\cdot d$$

$$\operatorname*{argmin}_{\{\mathbf{r}_{k}\}}\sum_{\rho}\Big{|}\Big{|}\mathcal{F}^{RS}(\rho)-\mathcal{F}^{V}(\rho)\Big{|}\Big{|}_{V}^{2}$$

$$\begin{array}{rl}& L_3=A_1\\ & L_5=A_1\cup A_2\\ & L_7=A_1\cup A_2\cup A_3\\ & L_{11}=A_1\cup A_2\cup A_3\cup B_1\end{array}$$

$$\{\rho \}=(\{\alpha \}\otimes \{{\alpha }^{\prime }{\}}_{l\le 2})\cup \{\beta \}$$

$$\begin{array}{rl}{E}_{xx}=& -\sum _{ij}^{occ}(ij|ij)\end{array}$$

$$E_C^{MP2}=-\sum _{ij}^{occ}\sum _{ab}^{virt}\frac{(ia|jb)[2(ia|jb)-(ib|ja)]}{{\epsilon }_a+{\epsilon }_b-{\epsilon }_i-{\epsilon }_j}$$

$$(ia|jb)\stackrel{RI-X}{=}\sum _{\beta {\beta }^{\prime }}{\mathcal{F}}_{\beta }^X(ia)V_{\beta {\beta }^{\prime }}{\mathcal{F}}_{{\beta }^{\prime }}^X(jb)$$

$$E_C^{RPA}=\frac{1}{2\pi }{\int }_0^{\infty }d\omega Tr[\ln (1-{\chi }_0(i\omega )\cdot v)+{\chi }_0(i\omega )\cdot v]$$

$${\chi }_0(\mathbf{r},{\mathbf{r}}^{\prime };i\omega )$$

$$\stackrel{RI}{\simeq }$$

$$\doteqdot$$

$$\big{[}\chi_{0}^{RI}(i\omega)\big{]}_{\beta\beta^{\prime}}=\sum_{kk^{\prime}}\;M_{\beta k}\;\chi_{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\omega)\;M_{\beta^{\prime}k^{\prime}}$$

$${\chi }_0({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};i\omega )=\sum _{\tau }{c}_{\tau }(\omega ){\chi }_0({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};i\tau )$$

$${\chi }_0({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};i\tau )=G^{<}({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};i\tau )G^{>}({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};-i\tau )$$

$$G^{<}({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};i\tau )$$

$$G^{>}({\mathbf{r}}_k,{\mathbf{r}}_{k^{\prime }};-i\tau )$$

$$\int_{0}^{\infty}d\omega\Big{[}\frac{1}{E-i\omega}+\frac{1}{E+i\omega}\Big{]}=\pi$$

$$\operatorname*{argmin}_{w_{k},z_{k}}\left[\int_{\ln(E_{min})}^{\ln(E_{max})}du\;\bigg{|}\bigg{|}\sum_{k}w_{k}\Big{[}\frac{1}{e^{u}-iz_{k}}+\frac{1}{e^{u}+iz_{k}}\Big{]}-\pi\bigg{|}\bigg{|}^{2}\right]$$

$$\operatorname*{argmin}_{w_{k}^{p},\tau_{p}}\left[\sum_{k}\int_{\ln(E_{min})}^{\ln(E_{max})}du\;\bigg{|}\bigg{|}\sum_{p}w_{k}^{p}e^{-\tau_{p}\,e^{u}}-\Big{[}\frac{1}{e^{u}-iz_{k}}+\frac{1}{e^{u}+iz_{k}}\Big{]}\bigg{|}\bigg{|}^{2}\right]$$

Full text

††preprint: AIP/DUCHEMIN

I Introduction

The resolution-of-identity (RI)Whitten (1973); Dunlap, Connolly, and Sabin (1979); Mintmire and Dunlap (1982); Van Alsenoy (1988); Feyereisen, Fitzgerald, and Komornicki (1993); Vahtras, Almlof, and Feyereisen (1993); Klopper and Samson (2002) stands as a central technique in quantum chemistry, relying on the expansion of $\phi_{n}\phi_{m}$ co-densities over an auxiliary atomic basis set $\{\beta\}$ that scales linearly with the number of atoms. Even though the auxiliary basis sets are typically three times larger than the corresponding atomic orbital (AO) basis sets supporting the $\phi_{n}$ Hartree-Fock or Kohn-Sham molecular orbitals, this represents a considerable saving both in terms of memory and number of operations to be performed, which comes at the price of a moderate accuracy loss.

While RI was introduced initially to facilitate the calculation of 2-electron 4-center $(mn|kl)$ Coulomb integrals, operators that finds an exact expression in the product space between occupied and virtual eigenstates can also be expressed more compactly using auxiliary bases. As such, the independent-electron density-density susceptibility $\chi_{0}({\bf r},{\bf r}^{\prime};\omega)$ is of central importance to the present study. In particular, the scaling of the related random phase approximation (RPA) approach, Bohm and Pines (1953); Gunnarsson and Lundqvist (1976); Langreth and Perdew (1977) a popular low-order perturbative approach to correlation energies beyond density functional approximations, can be reduced from $\mathcal{O}(N^{6})$ to $\mathcal{O}(N^{4})$. Eshuis, Yarkony, and Furche (2010)

The computational efficiency and accuracy of RI techniques strongly depend on the scheme adopted to build the appropriate coefficients expressing the $\phi_{n}\phi_{m}$ co-densities on the auxiliary basis. The original density-fitting (RI-SVS) approachDunlap, Connolly, and Sabin (1979); Van Alsenoy (1988); Feyereisen, Fitzgerald, and Komornicki (1993) expresses these coefficients as a direct overlap $\langle\phi_{m}\phi_{n}|\beta\rangle$, requiring the calculation of the sparse $\langle\alpha{\alpha}^{\prime}|\beta\rangle$ coefficients, with $\{\alpha\}$ being the AO basis set used to expand the molecular orbitals $\phi_{n}$. The now widely adopted Coulomb-fitting (RI-V) approximation Mintmire and Dunlap (1982); Vahtras, Almlof, and Feyereisen (1993) requires on the other hand the calculation of the denser 3-center Coulomb integrals $(\alpha{\alpha}^{\prime}|\beta)$, displaying much less sparsity than $\langle\alpha{\alpha}^{\prime}|\beta\rangle$ integrals due to the long-range nature of the Coulomb operator. The Coulomb-fitting formalism is known to be more accurate than the density-fitting scheme for auxiliary basis sets of similar sizes, Ren et al. (2012); Duchemin, Li, and Blase (2017) bringing to the standard issue of the trade-off between accuracy and computational/memory costs. The use of short-range or attenuated Coulomb operators, Jung et al. (2005); Reine et al. (2008); Sodt and Head-Gordon (2008); Ihrig et al. (2015) or corrective techniques such as “multipole-preserving” constraints to the density-fitting scheme,Van Alsenoy (1988); Duchemin, Li, and Blase (2017) allows one to tune the accuracy-to-cost ratio between these two standard RI approximations.

In conjunction with other powerful techniques, such as the Laplace transform, Almlöf (1991) exploiting the sparsity of the RI-SVS density-fitting $\langle\alpha{\alpha}^{\prime}|\beta\rangle$ coefficients in the limit of large systems was shown to allow cubic-scaling RPA calculations. Wilhelm et al. (2016) As a trade-off between accuracy and efficiency, a Coulomb-attenuated variation of the Coulomb-fitting (RI-V) RPA was recently explored to obtain a low-scaling formulation, Luenser, Schurkus, and Ochsenfeld (2017) exploiting further the decay properties in real-space of the Laplace-transformed ”pseudo” density matrices expressed in the AO basis. Häser (1993); Ayala and Scuseria (1999); Kállay (2015); Schurkus and Ochsenfeld (2016) The efficiency of this latter family of approaches depends on the electronic properties of the system of interest, and are different in nature from the sparsity associated with specific resolution-of-identity formalisms.

In the present study, we explore an alternative approach for reducing computational and memory loads by assessing on a large set of molecules the merits of a separable resolution-of-the-identity formalism relying on a density fitting scheme over compact sets of real-space points $\{{\bf r}_{k}\}$. Our approach preserves the use of standard Gaussian atomic orbitals and auxiliary basis sets for all-electron calculations, targeting the accuracy of the Coulomb-fitting (RI-V) formalism. The set-up of the fitting procedure scales cubically with system size. The accuracy of our approach is first validated by an extensive benchmark of the exchange and MP2 correlation energies over a large set of molecules. Combined with the Laplace transform technique, and following the so-called space-time approach for calculating the susceptibility operator,Rojas, Godby, and Needs (1995); Kaltak, Klimes̆, and Kresse (2014) the calculation of the RPA correlation energy within the present real-space quadrature approach is shown to scale cubically in terms of operations and quadratically in terms of memory, without invoking any sparsity or localization properties. The accuracy of the present real-space RI-RPA formalism is further shown to match that of the standard quartic-scaling Coulomb-fitting RI-RPA calculations for a large set of molecules including the oligoacenes, $C_{60}$ and a larger octapeptide angiotensin II molecule (146 atoms including 71 H atoms) proposed by Eshuis and coworkers in the early days of RI-RPA implementations. Eshuis, Yarkony, and Furche (2010)

II Theory

In this Section, we briefly outline the standard RI-V and RI-SVS approximations, introducing the notations used throughout the paper. We then discuss separable resolution-of-identities and present our specific implementation preserving the use of standard atomic and auxiliary basis sets. The present approach relies on weighted real-space $\delta({\bf r}-{\bf r}_{k})$-functions to express the density fitting coefficients, relating co-densities to auxiliary basis functions through real-space quadratures. The scheme to optimize the distribution of ${\bf r}_{k}$ and related weights is presented, and compared to other real-space quadrature formalisms. We demonstrate in particular that the computational cost associated with the setup of the present RI approach scales cubically with the system size. We show then how such a separable RI allows to obtain cubic-scaling RPA with low crossover with respect to standard quartic RI-RPA formalism when combined with the Laplace transform technique. We conclude this Section by presenting the technical details and parameters adopted in this study to perform the calculations illustrating the accuracy and scaling properties of the present approach.

II.1 Standard Resolution of the Identity

The resolution-of-identity (RI) approximationWhitten (1973); Dunlap, Connolly, and Sabin (1979); Mintmire and Dunlap (1982); Van Alsenoy (1988); Feyereisen, Fitzgerald, and Komornicki (1993); Vahtras, Almlof, and Feyereisen (1993); Klopper and Samson (2002) relies on the expansion of molecular orbitals co-densities $\phi\phi^{\prime}$ over an auxiliary basis set $\{\beta\}$, namely:

[TABLE]

The fit $\mathcal{F}$ is realized through an ensemble of measures $\{\mathcal{F}_{\beta}\}$, mapping the $\phi\phi^{\prime}$ product-space to the $\beta$ auxiliary subspace defined so as to scale linearly with the number of atoms. Typical examples of such procedures are the standard RI-V and RI-SVS fitting approaches that use respectively:

[TABLE]

where $V$ and $S$ represent respectively the Coulomb $(\beta|\beta^{\prime})$ and overlap $\langle\beta|\beta^{\prime}\rangle$ matrices associated with the auxiliary basis set, and $[X^{-1}]_{\beta\beta^{\prime}}$ denotes the $(\beta,\beta^{\prime})$ entry of the $X$ inverse matrix. To explicitly define our $\langle\cdot|\cdot\rangle$ and $(\cdot|\cdot)$ notations, we write:

[TABLE]

As shown in Ref. 13, both fitting techniques can be combined, preserving the Coulomb-fitting RI-V approach for low angular momentum auxiliary $\beta$ atomic orbitals. As emphasized here above, the number of $\langle\beta|\alpha\alpha^{\prime}\rangle$ overlap matrix elements in the RI-SVS approximation scales linearly with system size, offering a first strategy for reducing computational cost and memory thanks to sparsity. On the contrary, the number of $(\beta|\alpha\alpha^{\prime})$ Coulomb integrals scales quadratically so that sparsity, or sparse tensor algebra, is difficult to exploit within the more accurate Coulomb-fitting (RI-V) approach.

II.2 Separable RI

A separable expression for the resolution-of-the-identity can be obtained through a set of separable measures $\{\langle f_{k}|\}$ on the co-densities, namely:

[TABLE]

where the coefficients of $M$ have yet to be defined. Though it is not the only option, a trivial way to obtain separable measures is to work with $\delta(\mathbf{r}-\mathbf{r}_{k})$ distributions centered on a set of $N_{k}$ real-space locations $\{\mathbf{r}_{k}\}$. Working with real-space (RS) measures, the $\mathcal{F}^{RS}$ density fitting procedure takes then the simple form

[TABLE]

The clear advantage of the separability is that the two molecular orbitals $\phi\phi^{\prime}$, originally entangled in e.g. the $\mathcal{F}_{\beta}^{V}(\phi\phi^{\prime})$ RI-V fitting coefficients through the $(\beta|\phi\phi^{\prime})$ Coulomb integrals, are now disentangled. This will prove crucial in the calculation of linear-response or perturbation theory related quantities where summations over occupied/virtual pairs have to be performed as shown below in the case of the calculation of the RPA correlation energy. We emphasize however that, while relying on discrete values of the molecular orbitals in real-space, the present approach remains a resolution-of-the-identity in the sens that physical continuous quantities such as the co-densities, linear-response operators (e.g. suceptibility), etc. are defined everywhere in space in terms of the $\beta$ auxiliary basis functions.

Amongst existing formalisms adopting real-space quadrature strategies, several studies focused directly on 2-electron Coulomb integrals. The chain-of-sphere (COSX) semi-numerical approach to exchange integrals,Neese et al. (2009) building on Friesner’s pioneering pseudo-spectral approach,Friesner (1991) develops only one of the two co-densities forming 2-electron Coulomb integrals over a real-space grid. Alternatively, the Least Square Tensor Hypercontraction (LS-THC) formalism Parrish et al. (2012) fully develops the 2-electron integrals as a quadrature over real-space grid points, with an $\mathcal{O}(N^{4})$ computational complexity associated with the establishment of the quadrature.

On the other hand, the Interpolative Separable Density Fitting (ISDF) from Lu and coworkers provides a $\mathcal{O}(N^{3}\log(N))$ separable fit tensor by using Fourier transform and random projection techniques to select the $\mathbf{r}_{k}$ points and define the corresponding auxiliary densities, with a proof-of-concept application to a simple model system.Lu and Ying (2015, 2016) The present work adopts as well a separable form for the fit tensor (Eq. 5) through separable measures along real-space positions, but differs in the way the real-space locations $\mathbf{r}_{k}$ and their associated weights $[M]_{\beta k}$ are constructed, leading to a $\mathcal{O}(N^{3})$ quadrature determination process. In addition, we preserve the auxiliary atomic basis sets in use in quantum chemistry. The accuracy and efficiency of our approach is further benchmarked on a large number of molecular systems.

II.3 Quadrature weight determination

Assuming that the $\mathbf{r}_{k}$ locations are known (see below), the determination of the $[M]_{\beta k}$ coefficients in Eq. 5 are obtained from the following minimization condition:

[TABLE]

Namely, we aim to equate the fitting functions of the present RI-RS approach with that of the Coulomb fitting RI-V scheme. The accuracy of the present RI-RS approach is thus targeting that of the RI-V approach. The set of test co-densities $\{\rho\}$ typically spans the $\{\alpha\}\otimes\{\alpha\}$ product-space, even though it can be adjusted depending on the problem being addressed. Interestingly, the present fitting scheme allows to recover at reduced cost the otherwise $\mathcal{O}(N^{4})$ LS-THC factorization (see Supporting Information SI sup ). The solution of Eq. 6 can indeed be achieved with an advantageous $\mathcal{O}(N^{3})$ computational complexity, as demonstrated now.

In order to detail our fitting procedure, we adopt the following matrix notations: $[D]_{k\rho}=\rho(r_{k})$, $[F]_{\beta\rho}=\mathcal{F}_{\beta}^{V}(\rho)$. Due to the localization properties of the atomic orbitals, the number of atomic orbital products scales linearly with system size, and thus the number $N_{\rho}$ of test codensities in the test set can be considered $\propto N_{\alpha}$. As a result, the matrices D ($N_{k}\times N_{\rho}$), F ($N_{\beta}\times N_{\rho}$) as well as the matrix M ($N_{\beta}\times N_{k}$) of Eq. 5 are all $\mathcal{O}(N^{2})$ tensors. The fit equation Eq. 6 can then be formulated as:

[TABLE]

which leads to the standard least-square estimator:

[TABLE]

involving only matrix multiplications and inversions. Computation of $(D\cdot D^{\dagger})^{-1}$ could prove problematic if done explicitly. The term $D\cdot D^{\dagger}$ is positive, but has not guarantee to be definite. On the other hand, due to the large number $N_{\rho}$ of test co-densities, application of the standard SVD technique to extract the pseudo-inverse based estimator leads to rather significant prefactors to the otherwise $\mathcal{O}(N^{3})$ pseudo-inverse procedure. We take a side approach by combining simple balancing and Tikhonov $L_{2}$ regularization. Tikhonov, Goncharsky, and Yagola (1995) We first balance the problem by normalizing the rows of $D$, writing $\widetilde{D}=d\cdot D$ where $d$ is a diagonal matrix and the diagonal terms $[\widetilde{D}\cdot\widetilde{D}^{\dagger}]_{kk}=1$. The pseudo inverse is then calculated as:

[TABLE]

where the $L_{2}$ regularization parameter $\epsilon$ is adjusted to a small value to maintain definiteness of the problem and ensure numerical stability of the inverse. We identified the value $\epsilon=4\times 10^{-7}$ as a reasonable parameter for double precision arithmetic and kept this value for all the results presented in the current work. The resulting final least square estimator is thus:

[TABLE]

which can be computed efficiently through standard numerical inversion techniques. We emphasize that while computing the $[M]_{\beta k}$ optimal coefficients, there is no need for keeping the 3-center Coulomb integrals or the associated $\mathcal{F}^{V}_{\beta}(\alpha\alpha^{\prime})$ coefficients: these can be computed once on-the-fly and discarded immediately, avoiding thus any extra memory consumption. In other terms, we never store explicitly the $F$ and $D$ matrices but only their $F\cdot D^{\dagger}$ and $D\cdot D^{\dagger}$ ($N_{k}\times N_{k}$) resulting products.

II.4 Real-space grids generation

In the present approach, the optimized $\{{\bf r}_{k}\}$ sets are generated for isolated atoms, once for every chemical species and their associated atomic basis sets. These atomic grids are then duplicated according to the molecule geometry to form the system-specific quadrature points. With Eq. 9 defining the optimal $M$ for a given $\{{\bf r}_{k}\}$ set, locations are adjusted so as to minimize the fit error for a single atom test co-densities, using the Coulomb metric:

[TABLE]

For the sake of keeping the optimization process relatively simple, we structured the $\{{\bf r}_{k}\}$ set over four different shells, each one replicated with a different number of radii. The only parameters of the optimization problem are thus the number of radii and their length. The four base shells taken here and denoted $A_{1}$, $A_{2}$, $A_{3}$ and $B_{1}$ are subsets of the Lebedev quadrature gridsLebedev (1975) (denoted here $L_{i}$ for the Lebedev grid of order i) in the sense that:

[TABLE]

The determination of the number of radii associated with each shell has been done through experimentation until satisfying configurations were obtained. Base shells are provided in the SI Tables S4-7, and the resulting atomic quadrature grids for atoms H, C, N and O in the SI Tables S8-11.sup Optimizing freely all $\mathbf{r}_{k}$ locations and not only the radii may leads to significant improvement with respect to the grid size / accuracy ratio. Providing such grids falls however outside the scope of the present work.

II.5 Technical details

Benchmark Hartree-Fock and MP2 calculations were performed on a standard set of 28 medium size organic molecules containing unsaturated aliphatic and aromatic hydrocarbons or heterocycles, aldehydes, ketones, amides and nucleobases. Such a test set was originally proposed by Thiel and coworkers Schreiber et al. (2008) for reference optical excitations calculations within e.g. coupled cluster, Schreiber et al. (2008); Silva-Junior et al. (2010) TD-DFT Jacquemin et al. (2009) and more recently Bethe-Salpeter Jacquemin, Duchemin, and Blase (2015); Bruneval, Hamed, and Neaton (2015); Krause and Klopper (2017) formalisms. We adopt the MP2/6-31Gd geometries supplied in Ref. 35.

The assessment of the scaling properties of the present real-space quadrature RPA implementation is further performed on the oligoacenes family from benzene to hexacene using the B3LYP cc-pVTZ geometries available in Ref. 41, complemented by the decacene, recently observed, Krüger et al. (2017) and the (hypothetical) octacene, both relaxed at the B3LYP/6-31Gd level. Finally, we consider the $C_{60}$ fullerene (B3LYP/6-311Gd geometry provided in the SI) and the octapeptide angiotensin II molecule originally proposed by Eshuis and coworkers. Eshuis, Yarkony, and Furche (2010)

All calculations are performed with input molecular orbitals generated at the (spherical) cc-pVTZ Dunning (1989) Hartree-Fock level using the NWChem package.Valiev et al. (2010) The corresponding (cartesian) cc-pVTZ-RI auxiliary basis Weigend, Köhn, and Hättig (2002) was adopted in all resolution-of-identity (RI) approaches (RI-SVS, RI-V and real-space quadrature RI-RS). For sake of comparison, Hartree-Fock exchange and MP2 correlation energies were calculated exactly, namely without any RI approximation, using the NWChem package as well. All calculations are performed without any frozen-core approximation.

The set $\{\rho\}$ of test co-densities can be adjusted depending on the needs, for example to match a specific subset of the wave functions co-densities. In the rest of this work, we adopt the following settings:

[TABLE]

for both the single atom $\{\mathbf{r}_{k}\}$ sets problem and the full system optimization of $M$ coefficients. Limiting the second $\{\alpha^{\prime}\}$ atomic orbital (AO) basis set to $s$, $p$ and $d$ orbitals allowed to speed up the computation while no significant change of accuracy was observed. In the minimization process, the weight on the $s$ and $p$ AO orbitals have also been stressed with respectively factors 4 and 2, so as to increase focus on the low-order multipole charge and dipole component of the co-densities. Inclusion of the $\{\beta\}$ auxiliary orbitals within the test set slightly improves regularity of the errors.

We adopt real-space $\{{\bf r}_{k}\}$ sets that contain typically 320 points per C, N and O atom, and 180 for hydrogen. This size corresponds to about 3 times the size of the corresponding cc-pVTZ-RI auxiliary basis set. In the present study, we do not seek to look for minimal grid sizes, showing here below that excellent accuracy and small crossover between RI-RS and RI-V can already be obtained with such parameters. Details about the optimized real-space $\{{\bf r}_{k}\}$ sets, optimized for the cc-pVTZ and cc-pVTZ-RI Gaussian basis sets following Eq. 6 and 10, are provided in the SI. sup

The Laplace transform (LT) RPA correlation energy calculations are based on time and frequency grids described in the Appendix together with convergence tests for the benzene correlation energy, other molecules being reported in the SI.sup The present RI calculations, including the standard RI-V, RI-SVS and the newly developed real-space RI-RS, with and without Laplace Transform, are performed with a specific pilot code building on the Coulomb integral libraries implemented in the Fiesta code. Jacquemin, Duchemin, and Blase (2015); Duchemin, Li, and Blase (2017); Li et al. (2016)

III Results

III.1 Assessing the accuracy of the optimized real-space grid : exact exchange and MP2 correlation energies

As a first accuracy test of the present real-space RI implementation, namely to assess the quality of the co-density fits, we calculate both the exact exchange energy:

[TABLE]

written here for a spin compensated system, and the Møller-Plesset (MP2) correlation energy:

[TABLE]

The exchange and MP2 energies are calculated using the RI expressions of the 2-electron integrals, namely e.g.:

[TABLE]

where X=SVS, V or RS depending on the selected scheme. Since we want to specifically address the accuracy of the fitting technique, we avoid at that stage using extra approximations such as Laplace transform techniques.

The results, namely the error as compared to exact calculations for the 28 Thiel’s set molecules, are provided in Fig. 1 for exchange energy calculations. Clearly, the $\{{\bf r}_{k}\}$ sets adopted in this study provide errors that are of the same magnitude as the targeted RI-V approximation, and much smaller than those obtained with the RI-SVS scheme. We recall that our real-space quadrature was optimized to reproduce the coefficients $\mathcal{F}_{\beta}^{V}$ of the RI-V formalism (see Eq. 6) so that it should not be expected that the RI-RS approach yields errors smaller than the RI-V.

The RI-RS results are provided with an ”error bar” that represents the variance of the error distribution when the molecules are rotated. Contrary to the $\{\beta\}$ Gaussian auxiliary basis, the atomic $\{{\bf r}_{k}\}$ set is not rotationally invariant. Clearly, however, such a variance remains marginal as compared e.g. to the differential of errors between the standard Coulomb-fitting (RI-V) and density fitting (RI-SVS) schemes.

We now turn to the MP2 correlation energies (Fig. 2). Again, we observe that the RI-RS quadrature does not degrade significantly the targeted RI-V MP2 correlation energies, with an error remaining lower than a few tens of $\mu$Hartree per atom. Similarly, the variance of the error remains small. We provide in the Inset of Fig. 2 the actual distribution of errors obtained for a thousand independent random orientations of the benzene molecule, reporting the related variance error bar. We can conclude from the present set of benchmark calculations that the real-space representation (RI-RS), with the adopted $\{{\bf r}_{k}\}$ distributions size, reproduce accurately the RI-V Coulomb integrals involving co-densities (products $\phi_{i}\phi_{j}$ in the exact exchange expression) and transition densities (products $\phi_{i}\phi_{a}$ in the MP2 energy formula) at the core of all explicitly correlated perturbative techniques.

III.2 Laplace transformed RPA

We now turn to the central application of the present study, namely the calculation of the correlation energy within the random phase approximation (RPA). Pines and Bohm (1952); Nozières and Pines (1958); Langreth and Perdew (1977) We show in particular that the real-space quadrature RI-RS, combined with the standard Laplace transform (LT) approach, Almlöf (1991) allows reducing the scaling with system size down to $\mathcal{O}(N^{3})$, instead of the $\mathcal{O}(N^{4})$ scaling of standard RI-RPA implementations,Eshuis, Yarkony, and Furche (2010) without invoking any localization nor sparsity arguments.

Following seminal papers,Furche and Voorhis (2005); Fuchs et al. (2005); Furche (2008) we start with the adiabatic-connection fluctuation-dissipation theorem (ACFDT) formula to the RPA correlation energy :

[TABLE]

where $v$ is the bare Coulomb operator and $\chi_{0}(i\omega)$ the independent-electron density-density susceptibility at imaginary frequency, that is for closed-shell systems:

[TABLE]

with $(j,a)$ indexing (occupied/virtual) molecular eigenstates. The construction of the $\chi_{0}^{RI}(i\omega)$ matrix according to Eq. 17 clearly scales as $\mathcal{O}(N^{4})$.

To discuss such scaling properties, we compare first in Fig. 3 the total computing time for calculating the RI-RPA correlation energy within the standard Coulomb-fitting approach (RI-V) and the novel RI-RS formalism, using the acene family from benzene to decacene as a test set. Calculations are performed using the cc-pVTZ AO and cc-pVTZ-RI auxiliary basis sets, together with a 12-point quadrature rule for the imaginary frequency axis integration (see Appendix). The corresponding correlation energies are provided in the Appendix for benzene, and in the SI for other acenes, demonstrating again the accuracy of the real-space approach as compared to the targeted RI-V approach.

We observe that the RI-RS scheme walltime becomes smaller than that of our RI-V implementation for acenes larger than naphthalene. We emphasize however that at that stage, both RI-V and RI-RS techniques offer the very same $\mathcal{O}(N^{4})$ scaling, differing only by the expression of the $\mathcal{F}_{\beta}(\phi_{j}\phi_{a})$ coefficients. This crossover is related to the effort coming from the $\mathcal{O}(N^{4})$ computation of the full $\mathcal{F}_{\beta}(\phi_{j}\phi_{a})$ coefficients set. Without any assumption on the sparsity, the RI-RS scaling relies only on dense algebra techniques and can thus be implemented efficiently.

In order to reduce this scaling, one thus needs to avoid explicit calculation of the full $\mathcal{F}_{\beta}(\phi_{j}\phi_{a})$ coefficients set. We achieve this by evaluating the independent-electron susceptibility directly in the real-space representation before transforming it back to the normal auxiliary basis representation:

[TABLE]

The second step consists in applying the well known Laplace transform (LT) technique Almlöf (1991) so as to first compute $\chi_{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\tau)$ in the time domain where its expression is separable Rojas, Godby, and Needs (1995); Kaltak, Klimes̆, and Kresse (2014) and transform it back to the frequency domain, using quadrature rules to form $\chi_{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\omega)$ (Eq. 20). Such scheme allows us to work with factorized expression of $\chi_{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\tau)$ (Eq. 21):

[TABLE]

introducing the propagators of the occupied states and of the unoccupied states, respectively:

[TABLE]

with $\tau>0$ and the zero of occupied/virtual electronic energy levels taken at the Fermi level. As a result of the decoupling of occupied and virtual states, the $G^{<}$ and $G^{>}$ propagators can be obtained with $\mathcal{O}(N^{3})$ operations. Further, the entire \big{[}\chi_{0}^{RI}(i\omega)\big{]}_{\beta\beta^{\prime}} matrix stems from a combination of Hadamard product (Eq. 21) and standard matrix operations (Eq. 19-20), yielding an overall $\mathcal{O}(N^{3})$ process.

We can now report in Fig. 3 the full calculation walltime associated with the RI-RS+LT approach. Laplace transform quadratures for each of the 12 imaginary frequencies are performed with a grid of 18 imaginary times, yielding fully converged RPA energies (see Appendix). With such running parameters, the RI-RS+LT approach becomes more efficient than the standard RI-V formalism for systems larger than anthracene, outperforming the RI-RS approach for molecules larger than pentacene.

To better assess the scaling properties associated with the resolution of Eq. 15 within standard and Laplace Transformed RI-RPA approaches, we single out the corresponding computation walltimes in Fig. 4. We assume in particular that fitted co-densities coefficients are already available in the case of standard approaches, so that standard RI-V and RI-RS computational loads are equivalent. To probe larger systems, the test set is extended with the $C_{60}$ fullerene and the original octapeptide angiotensin II molecule proposed by Eshuis and coworkers.Eshuis, Yarkony, and Furche (2010)

Without accounting for the fitting of co-densities calculation, the Laplace Transform RI-RS RPA algorithm supersedes the standard RI-RS (or RI-V) RPA approach for systems larger than hexacene. This delayed crossover (as compared to Fig. 3) can be imputed to both the Laplace transform overheads and the fact that with the ${\{{\bf r}_{k}\}}$ real space point set presently used, \big{[}\chi_{0}(i\omega)\big{]}_{kk^{\prime}} matrices are roughly $3\times 3$ times bigger than \big{[}\chi_{0}(i\omega)\big{]}_{\beta\beta^{\prime}} ones, leading to a $3^{3}$ prefactor in the linear algebra operations. This last fact demonstrates the importance of operating on small ${\{{\bf r}_{k}\}}$ quadrature sets. The additional cost coming from obtaining the fitting parameters $\mathcal{F}_{\beta}(\phi_{j}\phi_{a})$ only adds to the cost of the no-Laplace-Transform RI approaches, bringing the overall crossover at the level of a pentacene molecule, as exemplified in Fig. 4. We finally observe $\mathcal{O}(N^{3.5})$ and $\mathcal{O}(N^{2.5})$ scaling laws (see dotted lines), indicating that the expected asymptotic $\mathcal{O}(N^{4})$ and $\mathcal{O}(N^{3})$ behaviours are not yet reached for the tested molecule sizes. We conclude this section by mentioning that the total walltime for the (RI-RS+LT) calculation of the cc-pVTZ RPA correlation energy for $C_{60}$ takes less than 6500 secs on a single processor.

IV Discussion

The present implementation can be compared to the real-space-grid imaginary-time approach introduced originally in the framework of $GW$ calculations by Rojas and coworkers Rojas, Godby, and Needs (1995) or the recent real-space-grid imaginary-time RPA implementation by Kaltak and coworkers. Kaltak, Klimes̆, and Kresse (2014) In such studies, the real-space grid was obtained as the Fourier transform of the planewave basis used to expend the Bloch states in a pseudopotential or PAW framework for periodic systems. Alternatively, the work of Moussa Moussa (2014) demonstrated as well cubic scaling for a random-phase approximation with second-order screened exchange formalism, exploiting a real-space grid for both the primary and auxiliary basis sets, combined with nested low-rank approximations to energy denominators. The present RI-RS approach, while targeting all-electron atomic-basis calculations, preserves the use of the standard representation of molecular orbitals, related co-densities and response operators in terms of atomic orbitals and their associated auxiliary bases, adopting further standard Laplace transform techniques.

A central issue in real-space representations, in particular when performing all-electron calculations, concerns the size of the real-space grid that strongly influences the crossover with standard RI implementations and the memory requirements. This is all the more important in the present study since we aim to calculate and store intermediate non-local operators such as the ${\chi}^{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\tau)$ susceptibilities, and not only local functions such as the charge density or the DFT exchange-correlation potential and energy density. Since our real-space $\{{\bf r}_{k}\}$ distribution must serve in a quadrature reproducing co-densities involving molecular orbital products, one may expect that it should be as large as standard grids Gill, Johnson, and Pople (1993) used to represent the charge density in DFT codes. Taking as an example the Gaussian09 code, the default DFT grid involves about 7000 grid points per atom after pruning. This is consistent with the recommended ”Grid3” in the original paper by Trutler and Alrichs Treutler and Ahlrichs (1995) yielding 5980 pruned points for elements from Li to Ne and that serves as the default in Turbomole.

Such standard DFT grid sizes are much larger than the number of real-space points used in the present RS approach, roughly 180 for Hydrogen and 320 per non-H atom (C, N, O). The present $\{{\bf r}_{k}\}$ sets were optimized so that the RI-RS scheme faithfully reproduces the standard RI-V density fitting results. Each set is composed of a number of different shells of high symmetry points, each one associated with a set of different radii. This minimization process results in a non-uniform ${\{{\bf r}_{k}\}}$ distribution of real-space points. As emphasized above, we did not seek to explore here in great details the set size to accuracy ratio, for our goal was to demonstrate that one can find an accurate real-space representation yielding a crossover with standard techniques for small system sizes. Smaller optimized sets of real-space points may certainly be explored in the future.

Another advantage the RI-RS+LT RPA scheme lies in the related $\mathcal{O}(N^{2})$ memory footprint associated with the underlying matrix algebra. In the case of the $C_{60}$ molecule, the relevant sizes are 1800 spherical AOs, 6060 cartesian auxiliary orbitals, and 19140 real-space quadrature points. On a single processor run, the set up of the RI-RS (Eq. 6) peaks at about 8 Gb of memory, while leaving at exit a memory footprint below 1Gb for the $M_{{\beta}k}$ coefficients. Further, each $\chi_{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\tau)$ requires about 3 Gb. In regards of a parallelization scheme oriented towards CPU efficiency, one can benefit from storing all the ($n\tau$) $\chi_{0}({\bf r}_{k},{\bf r}_{k^{\prime}};i\tau)$ in memory.

As emphasized here above, the present cubic scaling in terms of floating point operations, and quadratic scaling in terms of memory load, was obtained without invoking localization nor sparsity considerations. In particular, our approach does not require the use of the density-fitting (RI-SVS) approach with its sparse 3-center overlap matrix tensor. However, another class of localization properties, based on the exponential decay in real-space of the one-body Green’s function in gaped systems,Schindlmayr (2000) can be easily combined with the present approach. These localization properties, that strongly depend on the electronic properties of the system of interest, are reminiscent of the low-scaling techniques based on local AO formulations in the treatment of MP2 Häser (1993); Ayala and Scuseria (1999) or RPA Kállay (2015); Schurkus and Ochsenfeld (2016); Luenser, Schurkus, and Ochsenfeld (2017) correlation energies. not Such additional considerations, together with the stochastic approach by Neuhauser and coworkers, Neuhauser, Rabani, and Baer (2013) may be easily combined and explored in the future to further reduce the memory and computing time.

V Conclusion

We have introduced a separable resolution-of-the-identity based on a real-space quadrature of co-densities. The efficiency of our approach relies on setting up an optimal and compact distribution of real-space points $\{{\bf r}_{k}\}$ allowing excellent accuracy, as exemplified in the case of the exact exchange energy and further the MP2 and RPA correlation energies, taking as a test case a large set of molecular systems. Our approach preserves the use of standard Gaussian atomic orbitals and related auxiliary basis sets for all-electron calculations, the real-space set of points being used as an intermediate representation. We demonstrate that such an approach leads to calculating RPA correlation energies with a cubic scaling in terms of operations, quadratic in memory, without invoking any localization nor sparsity considerations that may be combined in the future. The limited number of needed real-space points allows early crossovers with traditional Coulomb-fitting RI-RPA calculations for systems as small as naphthalene or anthracene (Fig. 3). The application of such a real-space separable RI to other explicitly correlated techniques, such as the $GW$ and Bethe-Salpeter equation (BSE) formalisms for calculating charged and neutral electronic excitations in molecular systems,Blase, Duchemin, and Jacquemin (2018) are currently under exploration.

Acknowledgements.

The authors thank Thierry Deutsch and Pierre-Francois Loos for their critical reading of the manuscript and Denis Jacquemin for running the benzene RPA calculations with the Turbomole code. This research used resources from the French GENCI supercomputing centers under project no. A0030910016.

Appendix A Imaginary frequency and time quadratures

We briefly outline the imaginary frequency and time quadrature strategies adopted for calculating the RPA correlation energy and for setting the Laplace transform. Relying on the poles structure of the susceptibility (Eq. 16), we seek for a quadrature that reproduces the contribution of all possible poles in the imaginary axis integral, using the exact relation:

[TABLE]

where the pole energy $E$ can vary from $E_{min}$, i.e. the electronic energy gap, to $E_{max}$ determined by the maximum transition energy between occupied and virtual energy levels. Even if an optimal solution of this problem could be in principle determined through the minimax fitting approach usually applied to Laplace transforms,Kaltak, Klimes̆, and Kresse (2014) we prefer here a somewhat numerically simpler least square formulation, namely we seek for $n{\omega}$ frequencies $z_{k}$ and weights $w_{k}$ corresponding to:

[TABLE]

Here, the log scale has been preferred so as to obtain uniform oscillations of the error as a function of the integration variable (see Fig. 5), since integrating directly over the energy would favor larger errors in the small energy range.Kaltak, Klimes̆, and Kresse (2014) Under such conditions, minimax and least square approaches should provide consistent results. The $z_{k}$ frequencies and their corresponding $w_{k}$ weights are used for the integration reported in Eq. 15. Convergence tests are provided below in Table 1 for benzene, other acenes being dealt with in the Supporting Information (SI).sup Concerning benzene, our reference RI-V RPA correlation energy falls within less than 8 meV/molecule (0.19 kcal/mol) as compared to that obtained with other codes Bruneval et al. (2016); TUR with identified differences in the treatment of preceding Hartree-Fock calculations (RI vs no-RI) and use of the auxiliary basis (cartesian vs spherical). Together with potential differences in the treatment of the energy integration, these variations result in very small energy differences (see SI Table S2 for details).sup

For the Laplace transform, the time grid $\{\tau_{p},p\in[1,n\tau]\}$ is set accordingly so as to minimize the Laplace transform errors for the specific quadrature $(z_{k},w_{k})$ points obtained above, namely:

[TABLE]

where the weights $\{w_{k}^{p}\}$ depend on the targeted $z_{k}$ frequency. Again, the log scale is used so as to allow a regular sampling of the error oscillations. We draw the reader attention on the fact that the least square approach is clearly at advantage over the minimax fitting approach here since defining an alternant is not possible in the case of simultaneous fits (one fit for each $z_{k}$ frequency). The results provided here below for the benzene RPA correlation energy demonstrates already $<10^{-4}$Ha convergence for grids containing as few as $n\omega=6$ imaginary frequencies and $n\tau=9$ times. For $n\omega=12$ and $n\tau=18$, the parameters used in this manuscript for actual calculations, the correlation energy is clearly converged.