Cubic scaling algorithms for RPA correlation using interpolative separable density fitting

TL;DR

This paper introduces a cubic scaling algorithm for RPA correlation energy calculation that leverages Cauchy's integral formula and interpolative separable density fitting to significantly reduce computational cost.

Contribution

The authors develop a novel cubic scaling algorithm for RPA correlation energy using Cauchy's integral formula and interpolative separable density fitting, improving efficiency over previous methods.

Findings

Achieves cubic computational scaling for RPA correlation energy calculations.

Uses a geometrically convergent quadrature rule for the integral.

Reduces computational cost via interpolative separable density fitting.

Abstract

We present a new cubic scaling algorithm for the calculation of the RPA correlation energy. Our scheme splits up the dependence between the occupied and virtual orbitals in $\chi^0$ by use of Cauchy's integral formula. This introduces an additional integral to be carried out, for which we provide a geometrically convergent quadrature rule. Our scheme also uses the newly developed Interpolative Separable Density Fitting algorithm to further reduce the computational cost in a way analogous to that of the Resolution of Identity method.

Equations

$$\rho (x)=\sum _{j=1}^{N_{\text{occ}}}|{\psi }_j(x){|}^2,$$

$$E=T_{\text{s}}+U_{\text{ext}}+U_{\text{H}}+E_{\text{xc}},$$

$$T_{\text{s}}$$

$$U_{\text{H}}$$

$$E_{\text{x}}^{\text{EX}}$$

$$E_{\text{c}}^{\text{RPA}}$$

$${\hat{\chi }}^0(x,y,i\omega )$$

$$-\underset{0}{\overset{\infty }{\int }}{e}^{{\epsilon }_{j}t}{e}^{-{\epsilon }_{k}t}{e}^{-i\omega t}\text{d}t=\frac{1}{{\epsilon }_j-{\epsilon }_k-i\omega },$$

$$\frac{1}{2\pi i}\underset{\mathcal{C}}{\int }\frac{1}{(\lambda -{\epsilon }_j+i\omega )(\lambda -{\epsilon }_k)}\text{d}\lambda =\frac{1}{{\epsilon }_j-{\epsilon }_k-i\omega },$$

$$\frac{f_j-f_k}{{\epsilon }_j-{\epsilon }_k-i\omega }$$

$$=\frac{1}{2\pi i}\underset{\mathcal{C}}{\int }\frac{1}{(\lambda -{\epsilon }_j+i\omega )(\lambda -{\epsilon }_k)}\text{d}\lambda .$$

$$⟨x\left|{\hat{\chi }}^0(i\omega )\right|y⟩$$

$$=\sum _j^{\text{occ}}\sum _k^{\text{vir}}\frac{{\psi }_j^{\ast }(x){\psi }_k(x){\psi }_k^{\ast }(y){\psi }_j(y)}{{\epsilon }_j-{\epsilon }_k-i\omega }+\text{c.c.}$$

$$=\frac{1}{2\pi i}\underset{\mathcal{C}}{\int }\left(\sum _j^{\text{occ}}\frac{{\psi }_j^{\ast }(x){\psi }_j(y)}{\lambda -{\epsilon }_j+i\omega }\right)\left(\sum _k^{\text{vir}}\frac{{\psi }_k(x){\psi }_k^{\ast }(y)}{\lambda -{\epsilon }_k}\right)\text{d}\lambda +\text{c.c.},$$

$${\psi }_j^{\ast }(x){\psi }_k(x)\approx \sum _{\mu =1}^{N_{\text{aux}}}{\psi }_j^{\ast }(x_{\mu }){\psi }_k(x_{\mu })P_{\mu }(x),$$

$${\psi }_j^{\ast }(x){\psi }_k(x)\approx \sum _{\mu =1}^{N_{\text{aux}}}{C}_{jk}^{\mu }{P}_{\mu }(x),$$

$${\hat{U}}_{\xi }^{\text{occ}}(x)$$

$${\hat{U}}_{\xi }^{\text{vir}}(x)$$

$$M_{st}(x)={\overline{\mathcal{U}}}_s^{\text{occ}}(x){\mathcal{U}}_t^{\text{vir}}(x),s=1,...,r_{\text{occ}},t=1,...,r_{\text{vir}},$$

$$QR=\left[\begin{array}[]{c}M\\ \text{conj}(M)\end{array}\right]E=\mathbf{M}E,$$

$$R_{N_{\text{aux}},N_{\text{aux}}}\ge \text{tol}\cdot R_{1,1}>R_{N_{\text{aux}+1},N_{\text{aux}+1}}.$$

$${\psi }_j^{\ast }(x){\psi }_k(x)\approx \sum _{\mu =1}^{N_{\text{aux}}}{\psi }_j^{\ast }(x_{\mu }){\psi }_k(x_{\mu })P_{\mu }^{\ast }(x).$$

$${\psi }_j^{\ast }(x){\psi }_k(x)\approx \sum _{\mu =1}^{N_{\text{aux}}}{\psi }_j^{\ast }(x_{\mu }){\psi }_k(x_{\mu })\text{Re}[P_{\mu }(x)].$$

$$⟨x\left|{\hat{\chi }}^0(i\omega )\right|y⟩$$

$$=\frac{1}{2\pi i}\underset{\mathcal{C}}{\int }\left(\sum _j^{\text{occ}}\frac{{\psi }_j^{\ast }(x){\psi }_j(y)}{\lambda -{\epsilon }_j+i\omega }\right)\left(\sum _k^{\text{vir}}\frac{{\psi }_k(x){\psi }_k^{\ast }(y)}{\lambda -{\epsilon }_k}\right)\text{d}\lambda$$

$$+\frac{1}{2\pi i}\underset{\mathcal{C}}{\int }\left(\sum _j^{\text{occ}}\frac{{\psi }_j(x){\psi }_j^{\ast }(y)}{\lambda -{\epsilon }_j-i\omega }\right)\left(\sum _k^{\text{vir}}\frac{{\psi }_k^{\ast }(x){\psi }_k(y)}{\lambda -{\epsilon }_k}\right)\text{d}\lambda$$

$$\approx \sum _{\mu \nu }\frac{1}{2\pi i}\underset{\mathcal{C}}{\int }[\left(\sum _j^{\text{occ}}\frac{{\overline{C}}_j^{\mu }{C}_j^{\nu }}{\lambda -{\epsilon }_j+i\omega }\right)\left(\sum _k^{\text{vir}}\frac{C_k^{\mu }{\overline{C}}_k^{\nu }}{\lambda -{\epsilon }_k}\right)$$

$$+\left(\sum _j^{\text{occ}}\frac{C_j^{\mu }{\overline{C}}_j^{\nu }}{\lambda -{\epsilon }_j-i\omega }\right)\left(\sum _k^{\text{vir}}\frac{{\overline{C}}_k^{\mu }{C}_k^{\nu }}{\lambda -{\epsilon }_k}\right)]\text{d}\lambda P_{\mu }(x)P_{\nu }(y)$$

$$\displaystyle=\sum_{\mu\nu}\frac{1}{2\pi i}\int\limits_{\mathcal{C}}\Bigg{[}\mathbf{J}_{\mu\nu}(\lambda,\omega)\mathbf{K}_{\mu\nu}(\lambda)+\overline{\mathbf{J}_{\mu\nu}(\overline{\lambda},\omega)}\overline{\mathbf{K}_{\mu\nu}(\overline{\lambda})}\Bigg{]}\,\text{d}\lambda\,\,P_{\mu}(x)P_{\nu}(y)$$

$$=\sum _{\mu \nu }{\chi }_{\mu \nu }^0(i\omega )P_{\mu }(x)P_{\nu }(y)$$

$$=⟨x\left|{\tilde{\chi }}^0(i\omega )\right|y⟩,$$

$${\mathbf{J}}_{\mu \nu }(\lambda ,\omega )$$

$${\mathbf{K}}_{\mu \nu }(\lambda )$$

$${\chi }_{\mu \nu }^0(i\omega )$$

$$-\Delta g=4\pi\Bigl{(}f-\int\limits_{[0,1]^{d}}f\Bigr{)}$$

$$(\hat{v}f)(x)=\underset{[0,1{]}^d}{\int }\hat{v}(x-y)g(y)\mathrm{d}y$$

$$\hat{v}(x-y)=\frac{1}{\pi }\sum _{\xi \ne 0,\xi \in {\mathbb{Z}}^d}\frac{1}{|\xi {|}^2}{e}^{2\pi i\xi \cdot (x-y)}.$$

$$|{\mathbf{B}}_{\mu }⟩=|{\mathbf{P}}_{\nu }⟩S_{\nu \mu }^{-1/2},$$

Full text

This work is partially supported by the National Science Foundation under grants DMS-1454939.

©  2017. This manuscript version is made available under the CC-BY-NC-ND 4.0 license

http://creativecommons.org/licenses/by-nc-nd/4.0/

1. Introduction

In Density Functional Theory (DFT) [10, 13], the ground state energy of a many-body quantum system is written as a functional of the density $\rho$. In the Kohn-Sham (KS) formalism of DFT [13], instead of considering the original interacting system of $N_{\text{occ}}$ electrons, we consider a system of $N_{\text{occ}}$ non-interacting electrons (the KS system) under a different external potential whose ground state density is identical to that of the interacting system. In this effective single-particle system, the ground state density is given by

[TABLE]

where $\{\psi_{j}\}$ are the Kohn-Sham orbitals, the eigenstates of the effective single-particle system. It is assumed throughout that the KS orbitals are ordered so that $\psi_{1}$ is the ground state of the KS system, $\psi_{2}$ is the first excited state, and so on. In KS-DFT, the ground state energy of a system with $N_{\text{occ}}$ interacting electrons can be written as

[TABLE]

where

[TABLE]

are, respectively, the kinetic energy of the effective single-particle system, the potential energy due to the external potential $V_{\text{ext}}$, and the so-called Hartree energy, which represents the classical contribution of the energy from the Coulomb interaction between electrons. The remaining term in (1.2), $E_{\text{xc}}$, is known as the exchange-correlation energy. It has no known simple form in terms of the density $\rho$ or the Kohn-Sham orbitals $\{\psi_{j}\}$ and therefore needs to be approximated. There are many ways [21] of approximating this functional. In this work, we consider one of the more accurate (and more computationally expensive) approximations, the Random Phase Approximation (RPA). In particular, we separate out the exchange and correlation parts: $E_{\text{xc}}=E_{\text{x}}+E_{\text{c}}$, and we use the exact exchange $E_{\text{x}}^{\text{EX}}$ for the exchange energy $E_{\text{x}}$, and the Random Phase Approximation to approximate the correlation energy $E_{\text{c}}$ [23].

[TABLE]

where

[TABLE]

and $\widehat{v}$ is the Coulomb kernel (in particular, we will consider the periodic Coulomb kernel, see Section 2.4), and $\operatorname*{tr}[A]=\int\limits\left\langle x\left|A\right|x\right\rangle\,\text{d}x$. We will only consider the zero temperature case. This means that, in the ground state, the first $N_{\text{occ}}$ KS orbitals are filled while the rest are empty. So, $f_{\ell}=1$ if $1\leq\ell\leq N_{\text{occ}}$ (the occupied orbitals), and $f_{\ell}=0$ if $\ell>N_{\text{occ}}$ (the virtual orbitals).

In practice, one needs a way to obtain the KS orbitals before the energy can be computed. This can be done via a self-consistent iteration. However, we do not consider this here. Instead, we only consider the calculation of the energy after the KS orbitals are known. In this sense, we are considering a perturbative, non-self-consistent calculation of the RPA correlation energy. That is, in a practical implementation, the KS orbitals could be calculated via a self-consistent iteration using a computationally less expensive, but also less accurate, approximation for the exchange-correlation energy functional (e.g., LDA, GGA). The orbitals which are output from that self-consistent iteration can then be used to compute a more accurate approximation to the true correlation energy by using them to calculate the RPA correlation energy. In this way, one obtains an approximation to the true correlation energy which is better than the less expensive method (LDA, GGA, etc.), but also does not require the self-consistent iteration to deal with the expense of RPA.

Of all the terms we have defined above, the RPA correlation energy $E_{\text{c}}^{\text{RPA}}$ is the most computationally expensive to calculate. The goal of this paper is to provide a cubic scaling algorithm for the computation of this term. Before we can effectively talk about scaling, we must first define some notation. In this work, we will use a spatial discretization with equally spaced grid points. We denote the total number of grid points by $n$. We denote the number of occupied orbitals, i.e., the number of electrons, by $N_{\text{occ}}$. Since there are infinitely many KS orbitals $\{\psi_{j}\}_{j=1}^{\infty}$ 111When the spatial discretization is fixed with $n$ grid points, the total number reduces to $n$, which is much larger than $N_{\text{occ}}$. and the orbitals corresponding to higher energies will tend to have smaller contributions to $\widehat{\chi}^{0}$, we choose to use only the $N_{\text{orb}}$ KS orbitals of lowest energy in the RPA calculation. This gives us $N_{\text{vir}}=N_{\text{orb}}-N_{\text{occ}}$ virtual orbitals. Since $n$, $N_{\text{occ}}$, $N_{\text{vir}}$, and $N_{\text{orb}}$ all grow linearly with the system size, we will sometimes refer to a general $N$ as a characterization of the system size.

Very recently, a few cubic scaling methods for calculating the RPA correlation energy have been presented in the literature. The general idea involved is to split up the $j$ and $k$ dependence in the computation of $\widehat{\chi}^{0}$ by introducing a new integral. The idea is easy to motivate. From (1.6), we can see that everything can be done in cubic scaling if we are able to construct the matrices $\widehat{v}$ and $\widehat{\chi}^{0}(i\omega)$ in cubic time. $\widehat{v}$ is not hard to construct, so we will focus on $\widehat{\chi}^{0}$. $\widehat{\chi}^{0}$ has $O(N^{2})$ entries, so each entry of $\widehat{\chi}^{0}$ must be calculated in $O(N)$ time. By inspection of (1.7), it is clear that the most natural computation will take $O(N^{2})$ due to the coupling of $j$ and $k$ in the denominator. But if we can decouple the $j$ and $k$ dependence, then we can sum over each index separately and calculate each entry of $\widehat{\chi}^{0}$ in $O(N)$. Two different integrals have been utilized for this purpose. The first is

[TABLE]

where $\varepsilon_{k}>\varepsilon_{j}$. Using this integral, one separates the dependence of $j$ and $k$ in (1.7) into a product of exponentials inside the integral. This leads to the Laplace transform cubic scaling methods. This idea was first applied to RPA calculations in [12] and [11], where a projector augmented wave (PAW) basis was utilized. The idea was later extended to atomic orbitals [24, 19, 26]. The other integral used to break up the $j$ and $k$ dependency is

[TABLE]

where $\mathcal{C}$ is a positively oriented closed contour that encloses $\varepsilon_{j}-i\omega$, but not $\varepsilon_{k}$. This idea was first used in the context of cubic scaling RPA in [20]. Our algorithm in this paper will also adopt this idea.

The most significant contribution of this paper is the reduction of the prefactor in front of the cubic term in the computational complexity. In order to motivate this second main idea of this paper, let us examine how the density fitting (also called resolution of identity) approximation lowers the computational cost in the quartic scaling method [7, 22]. The idea behind the approximation is that $\widehat{\chi}^{0}$ is nearly equal to a low rank matrix due to its structure. So, $\widehat{\chi}^{0}$ (and $\widehat{v}$) are formed into smaller sized matrices (by writing $\widehat{v}$ into a smaller auxiliary basis and $\widehat{\chi}^{0}$ into the dual basis) before the trace is taken. The smaller matrix sizes lower the computational cost, but it is still $O(N^{4})$ since the coupling of $j$ and $k$ is unaffected by the approximation.

After we split up $j$ and $k$ in the denominator of (1.7) using Cauchy’s integral formula, we wish to further reduce the computational cost by taking advantage of the “low rank” nature of $\widehat{\chi}^{0}$ using the same idea as density fitting (DF). However, the DF approximation cannot be used in our case for two reasons. The first is that the DF itself takes $O(N^{4})$ operations, which destroys the cubic scaling. The second problem is that $j$ and $k$ are coupled in the coefficients of the density fitting method. As long as $j$ and $k$ remain coupled, $\widehat{\chi}^{0}$ cannot be constructed in $O(N^{3})$. Solutions to both of these problems are provided by the interpolative separable density fitting (ISDF) method [18]. The ISDF is capable of computing a decomposition of $\psi_{j}^{*}\psi_{k}$ similar to that of DF except that the $j$ and $k$ dependence in the coefficients are separated. Additionally, the decomposition can be performed in $O(N^{3})$ due to the use of a random projection in the method. Our use of the ISDF also reduces the memory cost of our algorithm to $O(nN_{\text{aux}})$ compared to $O(N_{\text{aux}}^{3})$ for the traditional resolution of identity approach.

Let us also mention the recent work [16], where a related problem of phonon calculation is approached from the point of view of the Sternheimer equations to represent $\widehat{\chi}$ acting on functions. Normally, for phonon calculations, $O(N^{2})$ Sternheimer equations would need to be solved in order to compute $\widehat{\chi}^{0}\widehat{v}$. However, by use of interpolative separable density fitting and a polynomial interpolation, they reduce the number of Sternheimer equations to $O(N)$, which enables a cubic scaling algorithm.

The rest of the paper is organized as follows. In Section 2, we reformulate the expression (1.6) into a new, approximate form. This new expression, characterized by (2.27), is used as the basis for our cubic scaling algorithm. Section 3 begins with a summary of our algorithm followed by a detailed description of each step. In Section 4, we run numerical tests to compare the scaling of our algorithm against the quartic scaling resolution of identity algorithm.

2. Derivation of the method

In this section, we reformulate the RPA correlation energy (1.6) which will be used to construct our cubic scaling algorithm. We will consider a computational domain with periodic boundary condition: Without loss of generality, up to a rescaling, the computational domain is taken to be $[0,1]^{d}$, where $d$ is the spatial dimension. For simplicity, in this paper we just consider the $\Gamma$-point calculation with periodic Coulomb kernel, while we will leave Brillouin-zone sampling to future works. The interpolative separable density fitting has been extended to Bloch waves in [18].

2.1. Contour integral representation

The first key idea is to split up the dependence of $j$ and $k$ in the denominator of (1.7). This is accomplished through the use of Cauchy’s integral formula. For a given $\omega$, let $\mathcal{C}$ be a closed contour in the complex plane oriented in the clockwise direction which encloses $\varepsilon_{\ell}$ for all $\ell$ which are unoccupied, and does not enclose $\varepsilon_{\ell}\pm i\omega$ for any $\ell$ that is occupied. An example of such a contour is shown in Figure 1.

While in principle the contour can be chosen differently for different $\omega$, later in Section 3.2, we will make the restriction that $\mathcal{C}$ is the same for all $\omega$ for the purpose of reducing computational costs. If $j$ is occupied and $k$ is unoccupied, then using Cauchy’s integral formula, we may write the coefficient in (1.7) as

[TABLE]

This can then be used to obtain an expression for $\widehat{\chi}^{0}$ which can be computed in cubic time,

[TABLE]

where c.c. is the complex conjugate. We show in Lemma 3.1 that the contour integral can be discretized with a number of quadrature points which is logarithmic in $(\varepsilon_{N_{\text{orb}}}-\varepsilon_{N_{\text{occ}}})/(\varepsilon_{N_{\text{occ}+1}}-\varepsilon_{N_{\text{occ}}})$.

Note that the formula (2.2) already provides a cubic scaling method for calculating $\chi^{0}(i\omega)$. In particular, ignoring logarithmic factors, $\chi^{0}(i\omega)$ can be calculated with cost $O(N_{\text{orb}}n^{2})$. However, the number of grid points $n$ could be much larger than the number of orbitals, so to reduce the prefactor of the computational cost, we will write the problem into an auxiliary basis set instead of using the spatial grid points.

Moreover, reducing the rank of $\widehat{\chi}^{0}$ from $N_{\text{occ}}\cdot N_{\text{vir}}$ (its rank before spatial discretization) to the number of grid points $n$ (its rank after spatial discretization, assuming $n<N_{\text{occ}}\cdot N_{\text{vir}}$) is really just a limitation placed on the operator by the particular discretization of the problem, rather than something inherent to the operator itself. The intuitive idea for the expected approximate low rank of $\widehat{\chi}^{0}$ is that $\widehat{\chi}^{0}$ contains $O(N_{\text{orb}})$ information in its definition, and therefore its approximate rank should scale linearly with the number of orbitals used in the calculation, rather than the number of grid points. This motivates the use of ISDF, as recalled in the next subsection.

2.2. Interpolative Separable Density Fitting

We use the Interpolative Separable Density Fitting (ISDF) [17, 18] to further accelerate the computation. ISDF aims at a representation of the orbital pair functions as

[TABLE]

where the $\{x_{\mu}\}$ and $\{P_{\mu}\}$ are chosen by the ISDF algorithm, which we will recall below for completeness of the presentation. Here $N_{\text{aux}}$ denotes the number of auxiliary orbitals needed to represent the orbital pairs involved; it is empirically established that $N_{\text{aux}}$ depends linearly on $N_{\text{orb}}$ [17], which we will also further verify in our numerical examples. The representation (2.3) should be compared to the traditional density fitting which yields

[TABLE]

where $\{P_{\mu}(x)\}$ are inputs to the DF algorithm and the coefficients $C_{jk}^{\mu}$ are determined via least squares fitting (in the $L^{2}$ or Coulomb metric). In (2.3), the $\psi_{j}^{*}(x_{\mu})\psi_{k}(x_{\mu})$ factor is the coefficient for the basis function $P_{\mu}(x)$. The main difference is thus that the $j$ and $k$ dependence of the coefficients are cleanly separated in ISDF, but not in DF. This is important for achieving cubic scaling in our work. Furthermore, ISDF has some other advantages over DF: the time and memory cost of ISDF is cheaper, in particular, it requires only $O(nN_{\text{aux}})$ memory and cubic scaling computational cost. In addition, the auxiliary basis functions are determined by the algorithm and do not have to be specified by the user.

Let us now describe the ISDF algorithm. The essential idea of the ISDF algorithm is to select important grid points $\{x_{\mu}\}$ via a randomized column selection algorithm. For the application to RPA correlation energy, we only need orbital pair functions of the type $\psi_{j}^{*}(x)\psi_{k}(x)$ where one of $j$ or $k$ is occupied and the other is unoccupied (see (2.2)). We can use this to our advantage by making a slight modification to the algorithm in [18], which would otherwise give an approximation for all $N_{\text{orb}}^{2}$ orbital pair functions. A version of ISDF which only calculates approximations for the orbital pair functions we are interested in is presented in Algorithm 1.

Compared with the original ISDF algorithm presented in [18], one technical difference is Step 5 in Algorithm 1. The reason for introducing $\mathbf{M}$, instead of just using $M$ in the QRCP there, is that now we can take the basis functions $\{P_{\mu}\}$ to be real. This can be seen by switching $j$ and $k$ in (2.3) and taking the complex conjugate,

[TABLE]

Note that both (2.3) and (2.9) are valid only because we included the conjugate of $M$ in the QRCP. Now, we may average the two expressions to show that we may write $\psi_{j}^{*}(x)\psi_{k}(x)$ in terms of real basis functions,

[TABLE]

It turns out that taking the basis functions to be real considerably simplifies the expression for $\chi^{0}$ that we will obtain. This leads to reduced computational effort as well as simpler code. For these reasons, we will always assume that the auxiliary basis functions from the ISDF are real.

2.3. Representation of $\widehat{\chi}^{0}$ using interpolative separable density fitting

Now we can use the ISDF to reduce the computational cost of the cubic scaling method for the RPA correlation energy. First, we approximate the $\widehat{\chi}^{0}$ operator by an operator $\widetilde{\chi}^{0}$ by using the ISDF approximation. For simplicity of notation, we define $C_{j}^{\mu}=\psi_{j}(x_{\mu})$.

[TABLE]

where the last line defines notation of $\widetilde{\chi}^{0}$, and we have used the short hands

[TABLE]

We note that in (2.11), the separability of the ISDF coefficients into the $j$ and $k$ components is crucial. Without this separability (e.g., if a conventional DF was used) we would not be able to calculate $\widetilde{\chi}^{0}$ in cubic time since the sums over $j$ and $k$ would not decouple.

Before continuing, we state explicitly our notation related to $\widehat{\chi}^{0}$ for the sake of clarity:

• •

$\widehat{\chi}^{0}$ – the original operator.

• •

$\widetilde{\chi}^{0}$ – the approximation to $\widehat{\chi}^{0}$ that is obtained by applying the ISDF approximation.

• •

$\chi^{0}$ – defined by (2.15). In (2.26), we will show that it is $\widetilde{\chi}^{0}$ in the dual basis to the auxiliary basis functions.

• •

The argument $i\omega$ is often suppressed below for simplicity of notation.

2.4. RPA correlation energy with auxiliary basis functions

We now return our attention to (1.6), where $\widehat{v}$ is the periodic Coulomb kernel (see e.g., [14]), which is defined such that $g=\widehat{v}f$ solves the Poisson equation

[TABLE]

with periodic boundary conditions and such that $\int\limits_{[0,1]^{d}}g=0$, to fix the arbitrary constant. $\widehat{v}$ can be written as an integral operator,

[TABLE]

where the kernel function, understood as a function in $L^{2}([0,1]^{d})$, is given by

[TABLE]

We want to rewrite the expression (1.6) using the approximate basis $\left\{\left|P_{\mu}\right\rangle\right\}$. However, since we are considering a problem with periodic boundary conditions, it will be advantageous for us to instead consider the basis $\left\{\mathbb{1},\left\{\left|\mathbf{P}_{\mu}\right\rangle\right\}_{\mu=1}^{N_{\text{aux}}}\right\}$, where $\left|\mathbf{P}_{\mu}\right\rangle$ is the shift of $\left|P_{\mu}\right\rangle$ so that it has zero mean, and $\mathbb{1}$ is the constant function with norm 1. The reasons for this change are explained further in Section 3.1. Since we wish to work with an orthonormal set, we introduce the orthonormalized basis functions

[TABLE]

where $S_{\mu\nu}=\left\langle\mathbf{P}_{\mu}|\mathbf{P}_{\nu}\right\rangle$. Then we can take the trace with respect to the orthonormal set

[TABLE]

In the following, we consider $\widehat{\chi}^{0}$ and $\widehat{v}$ to be linear operators on an $n$-dimensional space (i.e., the discretization of the operators with respect to our spatial grid) for the purposes of justifying our steps.

[TABLE]

The first line is justified as follows. First, we note that $\widehat{v}$ is bounded on finite dimensional spaces. Therefore, for $\widetilde{\chi}^{0}$ close enough to $\widehat{\chi}^{0}$, we have $\left\lVert\widehat{\chi}^{0}\widehat{v}-\widetilde{\chi}^{0}\widehat{v}\right\rVert\ll 1$. Thus, assuming that $I-\widehat{\chi}^{0}\widehat{\nu}$ is invertible and $\ln(I-\widehat{\chi}^{0}\widehat{\nu})$ makes sense, the following linear approximation is justified

[TABLE]

where $A:B$ means $\sum_{i}\sum_{j}A_{ij}B_{ij}$, i.e., the sum of the entries of the entrywise product. Before continuing our derivation, we first give a series expression for $\ln(I-\widetilde{\chi}^{0}\widehat{v})$. To justify the expansion, we assume that the eigenvalues of $\widehat{\chi}^{0}\widehat{v}$ are contained in the left half complex plane. Then for $\widetilde{\chi}^{0}$ sufficiently close to $\widehat{\chi}^{0}$, there exists $c>1$ such that the following expansion holds.

[TABLE]

The first thing to note about (2.22) is that the nullspace of $\widehat{v}$ is contained in the nullspace of $\ln(I-\widetilde{\chi}^{0}\widehat{v})$. Since $\widehat{v}$ is the periodic Coulomb operator, the constant function $\left|\mathbb{1}\right\rangle$ is in its nullspace (as $\xi=0$ is excluded from the summation in (2.18)). Therefore, we may drop the final term in (2.20), since it is zero. We remark that to accelerate the convergence with respect to the computational domain, a more sophisticated treatment of the Coulomb singularity at $\xi=0$, rather than taking the periodic Coulomb kernel, is often used. For example, see the method developed in [8].

Next, we note that the infinite sum in (2.22) is absolutely convergent, so we can continue (2.20) by applying the above expansion and taking the trace through the sum.

[TABLE]

Next, we write the terms $\sum_{\beta}\left\langle\mathbf{B}_{\beta}\left|(\widetilde{\chi}^{0}\widehat{v})^{p}\right|\mathbf{B}_{\beta}\right\rangle$ into a more computationally efficient, but approximate, representation. For simplicity, we only give the derivation for $p=1$, but the others are similar.

[TABLE]

where

[TABLE]

is the dual basis to $\left\{\left|\mathbf{P}_{\mu}\right\rangle\right\}$. Before commenting on the significance of this new representation, let’s write $\left\langle\widetilde{\mathbf{P}}_{\mu}\left|\widetilde{\chi}^{0}\right|\widetilde{\mathbf{P}}_{\nu}\right\rangle$ into a simpler form.

[TABLE]

where $\chi^{0}_{\mu\nu}(i\omega)$ is defined in (2.15), and can therefore be computed in cubic time. Let us define $v_{\mu\nu}=\left\langle\mathbf{P}_{\mu}\left|\widehat{v}\right|\mathbf{P}_{\nu}\right\rangle$. Then the right hand side of (2.24) reads $\operatorname*{tr}[\chi^{0}(i\omega)v]$, where the right hand side is just the standard trace of the product of the two matrices, $\operatorname*{tr}[AB]=\sum_{\mu}\sum_{\nu}A_{\mu\nu}B_{\nu\mu}$. Plugging this into (2.23), we obtain our final desired approximation,

[TABLE]

3. Algorithm

In this section, we present the cubic scaling algorithm for the calculation of the RPA correlation energy. We present a brief overview in Algorithm 2 before going into the details of each step. The computational effort is stated to the right of each step.

Without using the ISDF, the algorithm would be essentially exactly the same, but with Step 1 removed and Step 3a replaced by (2.2). Except, in the computational costs, each $N_{\text{aux}}$ would be replaced by $n$. So clearly, if $N_{\text{aux}}$ is much less than $n$, then including the ISDF can substantially speed up the algorithm.

3.1. Step 2 – Computing the Coulomb matrix

We note that $v$ can be efficiently computed by noticing that

[TABLE]

where

[TABLE]

Therefore, $\phi_{\nu}$ solves the Poisson equation, $-\Delta\phi_{\nu}=4\pi\mathbf{P}_{\nu}$ with periodic boundary conditions. This equation can be efficiently solved using the fast Fourier transform. Therefore, the functions $\phi_{\nu}$ can be precalculated at a total cost of $O(N_{\text{aux}}n\log n)$. Then the quadrature for (3.1) is straightforward.

We have glossed over a couple details here. First, the Poisson equation with periodic boundary conditions is not solvable unless $\mathbf{P}_{\nu}$ has an average value of 0. This is of course true by our construction of $\mathbf{P}_{\nu}$, and this is the reason for the use of the $\{|\mathbf{P}_{\mu}\rangle\}$ basis rather than the $\{|P_{\mu}\rangle\}$ basis when calculating the trace in (2.20). The second detail we’ve skipped is that the solution $\phi_{\nu}$ is not unique as we can add any constant and get another solution. However, it turns out that adding a constant to $\phi_{\nu}$ does not change the integral in (3.1) since $\mathbf{P}_{\mu}(x)$ has mean 0. So, this problem is also avoided by the use of the $\{|\mathbf{P}_{\mu}\rangle\}$ basis rather than the $\{|P_{\mu}\rangle\}$ basis.

3.2. Step 3 – Quadrature rule for the contour integral

Before discussing our proposed quadrature rule, let us discuss the symmetry of (2.11). For notational purposes, define

[TABLE]

It is straightforward to show the following symmetry across the real line,

[TABLE]

We wish to calculate

[TABLE]

where $\mathcal{C}$ is oriented clockwise and encloses $\{\varepsilon_{k}\}_{k=N_{\text{occ}}+1}^{N_{\text{orb}}}$ while enclosing none of $\{\varepsilon_{j}\pm i\omega\}_{j=1}^{N_{\text{occ}}}$. An example of such a contour is given in Figure 1. In order to use the symmetry about the real axis, we will enforce our contour to be symmetric about the real axis. Define $\mathcal{C}_{\text{upper}}$ and $\mathcal{C}_{\text{lower}}$ to be the parts of the contour in the upper and lower half plane. Then using (3.5) and (3.6), we have

[TABLE]

We construct a quadrature rule for Step 3a by using similar ideas as in [9, 15]. For simplicity, let us assume that $\varepsilon_{N_{\text{occ}}}=0$. To account for the fact that it is not, we will just have to shift the resulting quadrature points by $\varepsilon_{N_{\text{occ}}}$. For both simplicity of notation and to follow [9] more closely, let us define $m=\varepsilon_{N_{\text{occ}}+1}-\varepsilon_{N_{\text{occ}}}$ to be the energy gap and $M=\varepsilon_{N_{\text{orb}}}-\varepsilon_{N_{\text{occ}}}$. Then we map the rectangle with vertices $\pm K$ and $\pm K+iK^{\prime}$, where $K$ and $K^{\prime}$ are the complete elliptic integrals [1]

[TABLE]

to the upper half plane via two consecutive transformations $t\mapsto u\mapsto z$ (see Figure 2).

[TABLE]

where $\text{sn}(t|k^{2})$ is the Jacobi elliptic function. The values of $K$ and $K^{\prime}$ can be found, e.g., via the ellipkkp function in the Schwarz-Christoffel Toolbox for MATLAB [5]. The Jacobi elliptic functions $\text{sn}(t)$, $\text{cn}(t)$, and $\text{dn}(t)$ are implemented in the ellipjc function in the same toolbox. One must be careful when using such functions as there are a few different common conventions for how to parameterize the Jacobi elliptic functions. We have been using the parameter $k$ while the Schwarz-Christoffel Toolbox uses the parameter $L=-\ln(k)/\pi$.

The idea for the quadrature rule is as follows. We ultimately wish to construct a quadrature rule in the $z$-plane, but since the function we are integrating is periodic and analytic, one can show via some numerical analysis [4, Section 4.6.5] that we can construct a geometrically convergent quadrature rule by using the trapezoid rule in the $t$-plane. Essentially, the idea is to map the $z$-plane to a periodic rectangle, apply the trapezoid rule in the periodic rectangle where it is known to converge geometrically, and then map the quadrature points in the $t$-plane back to the $z$-plane to obtain our desire quadrature rule. The trapezoid rule has the added bonus of having a nice nesting property of the quadrature points, so that we can create an adaptive quadrature rule.

The numerical analysis tells us that the rate of convergence will be greatest when the quadrature path in the $t$-plane is as far away as possible from all singularities of the function we are integrating. In order to find the singularities in the $t$-plane, let’s first examine them in the $z$-plane. Due to symmetry, we will only consider the contour and singularities in the upper half plane. In our case, the function we are integrating is given by (2.11). For a given $\omega\geq 0$, its singularities (in the upper half plane) are $\{\varepsilon_{j}+i\omega\}_{j=1}^{N_{\text{occ}}}$ and $\{\varepsilon_{k}\}_{k=N_{\text{occ}}+1}^{N_{\text{orb}}}$. We first notice that the singularities we wish to avoid depend on $\omega$. There is nothing inherently difficult about this, and we could construct a different quadrature rule for each $\omega$. However, in order to save on computation, we want the quadrature rule to remain the same for each $\omega$. This way, $\mathbf{K}(\lambda)$ does not need to be recalculated for each $\omega$. Therefore, when we construct our quadrature rule, we wish to avoid all such singularities for $\omega\geq 0$. Since we have assumed that $\varepsilon_{N_{\text{occ}}}=0$, this implies that $\{\varepsilon_{j}+i\omega\}_{\omega\geq 0,j=1,...,N_{\text{occ}}}$ is contained in the left half $z$-plane. Therefore, it is sufficient for us to say that we wish to avoid the entire left half $z$-plane.

However, when we construct the quadrature rule, we are concerned with the singularities in the $t$-plane. By inverting the above mappings, we can see in Figure 2 that the left half $z$-plane is mapped to the region bounded by the yellow, green, and purple curves in the $t$-plane. So, it is sufficient for us to avoid this region. Next, we note that the rest of the singularities in the $z$-plane are contained in the red line. This line is mapped to the bottom edge of the rectangle in the $t$-plane. Finally, we wish for our contour in the $z$-plane to encircle the red line in the clockwise direction. This is achieved by taking a contour in the $t$-plane which goes from the left side of the rectangle to the right side. It is now clear how to maximize the distance between the singularities and the contour in the $t$-plane. We must draw our contour in the $t$-plane halfway between the bottom of the rectangle and the lowest point of the purple curve. This is demonstrated by a black horizontal line with X’s in Figure 2(a). We apply the trapezoid rule on this line. The line can be mapped back to the $z$-plane to create a quadrature rule as shown in Figure 2(c).

Rather than now going into the rigorous details of the above argument, we will simply state the results. The details and proofs are deferred to the Appendix. First, the algorithm for Step 3a is summarized in Algorithm 3.

Next, we state the convergence rate of the proposed quadrature rule, whose proof may be found in the Appendix.

Lemma 3.1**.**

Let $N_{\lambda}$ denote the number of quadrature points used to discretize (3.14) via the trapezoid rule. Then, for any $M/m>1$, the error of the quadrature rule is

[TABLE]

Therefore, our quadrature rule for the contour integral converges geometrically in the number of quadrature points. Additionally, the number of points $N_{\lambda}$ needed for convergence to a specific error tolerance increases only logarithmically as $(\varepsilon_{N_{\text{orb}}}-\varepsilon_{N_{\text{occ}}})/(\varepsilon_{N_{\text{occ}+1}}-\varepsilon_{N_{\text{occ}}})\to\infty$.

Finally, we note that for Step 3b, rather than calculating the trace of the log, it is more efficient to calculate the log of the determinant,

[TABLE]

This expression is nice for practical computation since it avoids the necessity of calculating the matrix logarithm. Additionally, the determinant of a matrix may be calculated easily via an LU decomposition, for which there are readily available scalable codes.

3.3. Step 4 – Quadrature rule for the frequency integral

For the $\omega$ integral, we used the following Clenshaw–Curtis quadrature [3, Eq 3.2] on the semi-infinite interval $[0,\infty)$,

[TABLE]

where

[TABLE]

where $L$ is a parameter that must be chosen (we used $L=10$). The value of $L$ can affect the number of grid points needed for convergence, but this dependence is not very sensitive. A necessary condition for the geometric convergence of this method is that $f(\omega)=O(\omega^{3/2})$ as $\omega\to\infty$ [3]. This is guaranteed by the following lemma.

Lemma 3.2**.**

$\operatorname*{tr}\left[\ln(I-\chi^{0}(i\omega)v)+\chi^{0}(i\omega)v\right]=O(\omega^{-2})$* as $\omega\to\infty$.*

Proof.

It is straightfoward to show that $|\chi^{0}_{\mu\nu}(i\omega)|<c_{\mu\nu}\omega^{-1}$, where $c_{\mu\nu}$ is independent of $\omega$. This implies $\left\lVert\chi^{0}(i\omega)\right\rVert_{F}<C\omega^{-1}$. Then we have

[TABLE]

Therefore, for $\omega$ large enough, the eigenvalues of $\chi^{0}(i\omega)v$ are all less than $1$ in magnitude. This justifies the Taylor series expansion in the following,

[TABLE]

where the constant $C$ changes from line to line. ∎

The use of Clenshaw–Curtis allows a simple and fast converging adaptive quadrature rule since the points of the quadrature rule have a nice nesting property as seen by (3.17).

4. Numerical results

Our numerical results use the following as the test problem. Our two dimensional spatial grid is $10\times 10N_{\text{occ}}$ equally spaced points. First, we solve for the KS orbitals of the periodic system with Hamiltonian $H=T+V$, where $T=-\frac{1}{2}\Delta$ is the kinetic energy operator and $V$ is the external potential. The external potential consists of $N_{\text{occ}}$ Gaussian potential wells, the centers of which are randomly perturbed from the centers of their respective $10\times 10$ box of grid points. Then the eigenvectors of $H$ are used as the orbitals in the calculation of the RPA correlation energy.

4.1. Convergence with respect to number of orbitals

In these tests, we check the convergence of the RPA energy with respect to the number of orbitals used in the calculation. Figure 4(a) shows the results for a system with 4 electrons (and therefore a maximum of 400 orbitals). In Figure 4(b), we have scaled the entire system up by a factor of 8 (maximum of 3200 orbitals) and run the same test. We can see by comparing the figures that the results are essentially identical. Both the percentage of orbitals needed in the calculation for a particular error value and the number of auxiliary basis functions (as a percentage of the number of grid points) needed for a particular error level in the ISDF step are nearly the same in the two cases.

In general, one would want to work in a regime where $N_{\text{aux}}$ is as small as possible while still achieving sufficient accuracy. Note that if $N_{\text{aux}}\approx n$, then there is no point in using ISDF and one would be better off using (2.2) instead. Figure 4 implies that ISDF is worth doing in the $10^{-2}$ relative error range, where we can see from Figure 4, the ISDF yields an $N_{\text{aux}}$ significantly below $n$. However, this statement is highly dependent on the number of grid points. For example, in Figure 5, we see that the ISDF can be worthwhile all the way down to $10^{-4}$ relative error in the $n=1600$ case. The reduction of the basis size from $n$ to $N_{\text{aux}}$ (as a proportion of the number of grid points) is amplified as the number of grid points is increased with all other variables fixed. This is because, as will be discussed shortly, $N_{\text{aux}}$ depends on $N_{\text{orb}}$, not $n$.

In Figure 5, we fix an external potential and look at the behavior of the algorithm when the number of primal basis functions (grid points) is increased. Both tests are run with $N_{\text{occ}}=4$. One is run with $n=400$ grid points and the other with $n=1600$. Therefore, the maximal number of auxiliary basis functions are 400 and 1600, respectively. We make two observations about the figure. First, the number of auxiliary basis functions required depends only on the number of orbitals $N_{\text{orb}}$ used in the calculation, until saturation occurs. That is, in both the $n=400$ and $n=1600$ tests, $N_{\text{aux}}$ is essentially the same for $N_{\text{orb}}\leq 200$, at which point the number of auxiliary basis functions starts to max out at 400 in the smaller system. Second, we note that the error in the two tests is mostly identical for a given number of orbitals $N_{\text{orb}}$ used in the calculation. These two observations suggest that the ISDF procedure behaves precisely how one would hope as the number of primal basis functions (grid points) is increased. That is, (before the auxiliary basis functions max out) the number of auxiliary basis functions and the error in $E_{\text{c}}^{\text{RPA}}$ depends only on $N_{\text{orb}}$, and not on $n$. This last observation is, of course, only valid when $n$ is large enough and tol in the IDSF is small enough that the errors due to the spatial discretization and the ISDF approximation are negligible compared to the error induced by truncating the number of orbitals.

4.2. Cubic scaling

In this test, we show the cubic scaling behavior of the algorithm. The quartic scaling method using traditional density fitting is also plotted for comparison. The traditional density fitting requires that we input basis functions, so we use the basis functions $\{P_{\mu}\}_{\mu=1}^{N_{\text{aux}}}$ from the ISDF. The results from Figure 4 suggest that as we scale up the system size, we can choose the number of orbitals $N_{\text{orb}}$ to use in the calculation as a constant percentage of the number of grid points. So, we choose $N_{\text{orb}}=0.2n$. We scale the system size up to a maximum of $N_{\text{occ}}=160$. We can see in Figure 6 that the cubic scaling algorithm greatly outperforms the quartic scaling algorithm for large system sizes.

5. Conclusion

In this paper, we have presented a new cubic scaling algorithm for the computation of the RPA correlation energy. The key of the algorithm was to separate the dependence on $j$ and $k$ in the denominator of (1.7). This allows a natural cubic scaling method. However, in order to further reduce the computational cost, we employed the ISDF in analogy to how density fitting is used in the quartic algorithm. Another key idea to keep the computational cost down was to take advantage of the periodic and analytic nature of the function in the contour integral, which resulted in a geometrically convergent nested quadrature rule based on the simple trapezoid rule.

It is worth noting that the algorithm presented is highly parallelizable. Step 1, the ISDF, is composed of linear algebra routines which can be parallelized. It is clear that Steps 2 and 3a can be parallelized. Step 3b is parallelizable using the comment at the end of Section 3.2. There is no need to parallelize Step 4 as it is a simple one dimensional integral, but the Step 3 computations for each quadrature point $\omega_{m}$ could also be done in parallel.

Future directions include a parallel implementation of the algorithm, as well as implementation into scientific software. Another direction would be to look into the analytic properties of $\chi^{0}(i\omega)$. It would also be interesting to apply the ISDF to the Laplace transform method for cubic scaling RPA algorithms. We also plan to extend the algorithm presented in this paper to particle-particle RPA (ppRPA) [25].

Appendix A Derivation of Algorithm 3

In this Appendix, we first rigorously derive Algorithm 3. Then we conclude by proving Lemma 3.1.

Recall from Section 3.2 that we wish to find the lowest point of the purple curve in the $t$-plane. We do this now. First, we note that (3.12) is a Möbius transformation and therefore its inverse maps the imaginary line to a generalized circle in the $u$-plane. In particular, it maps the upper half imaginary line in the $z$-plane to the upper semicircle with radius $k^{-1}$ centered at the origin (the purple curve in Figure 2(b)). To map this semicircle back to the $t$-plane, we note the formula for the inverse of $\text{sn}(t)$, [2, Chapter 11.3]

[TABLE]

Then we can use basic calculus to minimize $\text{Im}\left[\text{sn}^{-1}(k^{-1}e^{i\theta})\right]$ over $0\leq\theta\leq\pi$.

[TABLE]

This expression is 0 if and only if the expression under the radical is nonpositive. Since the imaginary part of the expression under the radical must be 0, we require $\theta\in\{0,\pi/2,\pi\}$. [math] and $\pi$ correspond to the corners of the rectangle, so this means that $\theta=\pi/2$ must give us the minimum imaginary part of points along the purple curve. In conclusion, we choose our quadrature points in the $t$-rectangle with imaginary part given by

[TABLE]

This integral must be carried out numerically. However, it is simple and only needs to be done once at the beginning of the algorithm, so we just use the midpoint rule. We also note that we can easily remove any guess work here by proving a practically useful bound which can be obtained via the standard error analysis for the midpoint rule. First let

[TABLE]

Then computation shows

[TABLE]

By noting $0\leq k<1$ and $0\leq sk\leq 1$, we have

[TABLE]

Using the fact that the right hand side of (A.6) is decreasing on $[0,\infty)$,

[TABLE]

where the last line uses the fact that the preceding line is a strictly increasing function of $k^{-1}$. This estimate implies that a mesh size of $1/100$ guarantees an accuracy of $10^{-4}$. Considering the fact that if the value of this integral is off by a little it will only slightly change the convergence rate, this is sufficiently accurate.

We conclude this discussion of the quadrature rule with some brief analytic results, including the proof of Lemma 3.1. First, we note that in a realistic system, $M/m\gg 1$ which implies $k\approx 1$. This guarantees that using the midpoint rule to calculate $I$ will not require more than about 100 grid points. Next, we note that $I>K^{\prime}/4$. This can be seen by showing that the circle with radius $k^{-1/2}$ centered at the origin in the $u$-plane maps to the horizontal line with imaginary part $K^{\prime}/2$ in the $t$-plane. Before proving this statement, let’s see why this implies $I>K^{\prime}/4$. First, note that $1<k^{-1/2}<k^{-1}$. Therefore, the circle with radius $k^{-1/2}$ in the $u$-plane must map between the purple and red curves in the $t$-plane. This means that the purple curve cannot go any lower than $K^{\prime}/2$, which implies $I>K^{\prime}/4$. To show that the aforementioned circle maps to a line with constant imaginary part, it is enough to show

[TABLE]

is equal to 0 for all $0\leq\theta_{1}\leq\theta_{2}\leq\pi$. It is easily verified that (for $0<k<1$) the denominator on the right is always positive and the numerator is always a pure imaginary number. Therefore, the integrand is always purely imaginary, which proves the claim. Finally, the imaginary part of the line is $K^{\prime}/2$ since $\text{sn}^{-1}(ik^{-1/2})=iK^{\prime}/2$ [6, Table 22.5.2].

Now following [9] and using the fact that $I>K^{\prime}/4$, we have that for any $M/m>1$, the error of the quadrature rule is

[TABLE]