Complex-time shredded propagator method for large-scale $GW$ calculations

The $GW$ method is a many-body electronic structure technique capable of generating accurate quasiparticle properties for realistic systems spanning physics, chemistry, and materials science. Despite its power, $GW$ is not routinely applied to study large complex assemblies due to the method's high computational overhead and quartic scaling with particle number. Here, the $GW$ equations are recast, exactly, as Fourier-Laplace time integrals over complex time propagators. The propagators are then “shredded” via energy partitioning and the time integrals approximated in a controlled manner using generalized Gaussian quadrature(s) while discrete variable methods are employed to represent the required propagators in real space. The resulting cubic scaling $GW$ method has a sufficiently small prefactor to outperform standard quartic scaling methods on small systems ($⪆10$ atoms) and offers 2–3 order of magnitude improvement in large systems ($\approx 200–300$ atoms). It also represents a substantial improvement over other cubic methods tested for all system sizes studied. The approach can be applied to any theoretical framework containing large sums of terms with energy differences in the denominator.

Figure 1

An example of the proposed energy windowing approach with $N_{a_w}=N_{b_w}=2$ (a) For gapped systems, the energy ranges of $\left\{a_i\right\}$ and $\left\{b_j\right\}$ do not overlap. (b) For systems with overlapping energy ranges, energy window pairs arise both with energy crossings, red arrows, and without, blue arrows.

Figure 2

Numerical error vs computational savings for our cubic scaling formalism, CTSP-W, compared to the standard quartic $GW$ formulation for bulk Si and MgO modeled in a 16 atom supercell. The CTSP-W error decreases and computational work increases as the integration error is decreased (i.e., the number of quadrature points is increased). Computational work is measured by the ratio of operation count, Eq. (16), of the cubic method to the quartic method. (Top) Error in the macroscopic optical dielectric constant $[{\epsilon }_{\infty }\left(\mathrm{MgO}\right)=6.35, {\epsilon }_{\infty }\left(\mathrm{Si}\right)=64.85]$. (Bottom) Error in the COHSEX band gap $[E^{(\mathrm{gap},\mathrm{COHSEX})}\left(\mathrm{MgO}\right)=7.56$ eV, $E^{(\mathrm{gap},\mathrm{COHSEX})}\left(\mathrm{Si}\right)=1.92$ eV].

Figure 3

Error in the macroscopic RPA optical dielectric constant ${\epsilon }_{\infty }$ for the interpolation, the CTSP-W and the CTSP-1 methods with respect to the prediction of the standard quartic $\mathcal{O}\left(N^4\right)$ technique. The horizontal axis is the ratio of computational load of the cubic methods to the standard $\mathcal{O}\left(N^4\right)$ method as measured by operation count for a system of 16 Si atoms. Upper: Bulk Si data generated by using input fractional quadrature error ${\epsilon }^{\left(q\right)}\{0.001,0.01,0.1,0.2\}$ for interpolation; $\{0.001,0.01,0.1,0.3,0.5\}$ for CTSP-1; and $\{0.001,0.01,0.1,0.3,0.5,0.8\}$ for CTSP-W. Middle: same for bulk MgO with ${\epsilon }^{\left(q\right)}\{0.001,0.01,0.1\}$ for interpolation $\{0.001,0.01,0.1,0.3,0.7\}$ for CTSP-1 and $\{0.001,0.01,0.1,0.2,0.4\}$ for CTSP-W. One point on each method's curve is generated per input value ${\epsilon }^{\left(q\right)}$.

Figure 4

Error in the COHSEX approximation to the band gap ($\mathrm{\Gamma }-X$ gap for Si and $\mathrm{\Gamma }$ for MgO) for the interpolation, and the CTSP-W with respect to the prediction of the standard quartic $\mathcal{O}\left(N^4\right)$ technique as function of the ratio of operation count to the standard method. All data sets are computed for a supercell size of 16 atoms. The nomenclature and numerical tolerances are those of Fig. 3.

Figure 5

The relation between the input fractional quadrature error tolerance ${\epsilon }^{\left(q\right)}$ and the output error in the physical observable ${\epsilon }_{\infty }$ for the CTSP-W method applied to Si and MgO. The dotted line represents the situation where quadrature error equals to observable error.

Figure 6

Error in the $P_{0,0}^{q=(\frac{1}{2},\frac{1}{2},\frac{1}{2})}$ element of bulk Al using CTSP-W. The horizontal axis is the ratio of computational load of the cubic to that of the standard $\mathcal{O}\left(N^4\right)$ method for a supercell of eight Al atoms. A total 400 states were used, and the broadening parameter was set to 0.03 Ry (see Appendix pp8). The set of fractional quadrature errors ${\epsilon }^{\left(q\right)}$ employed to generate the curve is $\{0.01,0.1,0.3,0.5,0.7,0.8\}$.

Figure 7

Error in the bulk $G_0{W}_0$ band gap ($\mathrm{\Gamma }-X$ gap for Si) for the CTSP-W method as a function of the ratio of computational load to that of the “exact” quartic method (horizontal dashed line). All data were generated using a supercell of 16 Si atoms. The set of fractional integration errors employed to generate the curve is $\{0.001,0.01,0.05,0.1,0.5,1\}$. The computed band gap with the $O\left(N^4\right)$ method is $E^{(\text{gap},G_0{W}_0)}=1.37$ eV.

Figure 8

Compute time per operation for evaluation of the static $P$ of crystalline Si. Black squares indicate the $\mathcal{O}\left(N^4\right)$ method, and red circles and blue asterisks indicate the $\mathcal{O}\left(N^3\right)$ CTSP-W method with maximum fractional quadrature error (${\epsilon }^{\left(q\right)}$) of 1% and 10%. A standard workstation with the Linux operating system was employed—the application was not parallelized. The flatness of the curves for both the $\mathcal{O}\left(N^4\right)$ and $\mathcal{O}\left(N^3\right)$ methods indicate their computation time increase as $N^4$ and $N^3$: there is no additional hidden cost to using the $\mathcal{O}\left(N^3\right)$ method in an actual calculation. Computational details are given in Appendix pp8. Issues such as cache utilization, different for different system sizes in our untuned application, account for the factors of $<2\times$ deviation from a horizontal line.

Figure 9

Comparison of the performance, absolute band gap error versus computational load, for the $\mathcal{O}\left(N^4\right) GW$ terminator method of Ref. [19] as implemented in the yambo software, labeled yambo BG, and the $\mathcal{O}\left(N^3\right)$ CTSP-W method for a 16-atom Si cell with one $k$ point. The BG band gap is shown as red circles. The CTSP-W 2464 band gap (i.e., computed using 2464 bands) is given by the green stars, and the CTSP-W 1280 band gap (i.e., computed using 1280 bands) is given by the blue squares. See Fig. 3 for other details.

Figure 10

Comparison between the performance of the cubic scaling Foerster et al. [22] and the CTSP-W method. The $y$ axis is the error (%) in $P^{\left(\mathrm{model}\right)}$ (for a 16-atom Si supercell with 399 bands) and the $x$ axis is the logarithm of the ratio of computational work to that of the quartic method. The standard for the error is the quartic scaling result. For the Foerster method, the operation count is $N_{\omega }(N_c+N_v)$ where $N_{\omega }$ is the number of frequency grid points, and for the CTSP-W method the workload is defined in Eq. (16) where $N_r$ is set to one. The band gap $E^{\left(\mathrm{gap}\right)}$ is manually adjusted (i.e., using a “scissors” operation, to investigate its effect on performance. F-single and F-double mean the Foerster method using single and double windows.

Figure 11

Comparison of the performance of the cubic scaling minimax method of Liu et al. [23] and the CTSP-W method of this paper, for the computation of $P^{\left(\mathrm{model}\right)}$. The computational work ($x$ axis) for the minimax method is $N_{\omega }(N_c+N_v),$ where $N_{\omega }$ is the number of imaginary time grid points (see also the text). The standard quartic scaling method provides the baseline.

Figure 12

Gauss-Laguerre (GL) quadrature error in the integration of $e^{-x\tau }$ with 12 quadrature points as a function of ${ln}_{10}x$, solid blue curve [see Eq. (A1)]. The dashed black horizontal line shows that for $-0.75⪅{log}_{10}x⪅0.75$ equal error is generated for $x$ and $1/x$.

Figure 13

The fractional error, ${\epsilon }^{\left(q\right)}$, of Gauss-Lauguerre quadrature (for $x=1/E^{\left(\mathrm{bw}\right)}$) as a function of $\alpha$ and $N^{(\tau ,\mathrm{GL})}$. Here, $\alpha$ is defined to be $\sqrt{E^{\left(\mathrm{bw}\right)}/E^{\left(\mathrm{gap}\right)}}$. Each contour is labeled by ${\epsilon }^{\left(q\right)}$. For a fixed fractional error, $N^{(\tau ,\mathrm{GL})}$ is linear in $\alpha$.

Figure 14

Computational cost to compute the static $P$ using a $N_{v_w}=2\times N_{c_w}=2$ window scheme as a function of energy window size for a 16-atom Si system with 399 states (32 occupied and 367 unoccupied states) and one $k$ point (no $k$-point sampling and hence $q=0$ strictly). The position of the minimum does not occur at the equipartition point, (1,1), but rather at $(E_v^{\left(\mathrm{ratio}\right)},E_c^{\left(\mathrm{ratio}\right)})=(1.25,0.29)$.

Figure 15

Minimized computational cost $C^{\left(\mathrm{GL}\right)}$ to compute the static $P$ of a 16-atom Si system for window pairs $\{N_{v_w},N_{c_w}\}$ spanning $(N_{v_w}=9\times N_{c_w}=9)$ (81 total pairs). The total number of bands in the system was taken to be 399 (32 occupied and 367 unoccupied states) and one $k$ point (no $k$-point sampling and hence $q=0$ strictly). The position of the minimum is at the point, $(N_{v_w}=1,N_{c_w}=5)$.

Figure 16

(Left) Comparison of the two weight functions described in the text for $\gamma =1$. The blue dashed curve is the exponential weight $exp(-\tau )$ associated with Lorentzian broadening; the solid red curve is the new weight function associated with Eq. (C3). (Right) Fourier transforms of the weight functions. The transform of the exponential weight $e^{-\tau }$ (dashed blue) is $x/(1+x^2)$ while the transform of the weight $h\left(\tau \right)=exp(-\tau -{\tau }^2/2)$ is given by Eq. (C3) (solid red). For comparison, the target function $1/x$ is shown as well (short dashed green). Equation (C3) is smoother for small $x$ and approaches $1/x$ more rapidly at large $x$ than $x/(1+x^2)$.

Figure 17

Hermite-Gauss-Laguerre quadrature error as function of $x$ and $N^{(\tau ,\mathrm{HGL})}$ [see Eq. (D1)]. (Top) The fractional quadrature error (${\epsilon }^{\left(q\right)}$) is indicated as blue dots along with the fit function, ${\epsilon }_{\mathrm{fit}}^{\left(q\right)}$ [see Eq. (D2)]. (Bottom) The contour lines of ${\epsilon }_{\mathrm{fit}}^{\left(q\right)}$ are shown. Each contour line can be represented with high fidelity using only quadratic function of $x$ [see Eq. (D3)].

Figure 18

The optimized cost function to compute the model dynamic ${\mathrm{\Sigma }}^{\left(\mathrm{model}\right)}\left(\omega \right)$ of the text as a function of the number of energy windows for bulk silicon with $\omega =E_{v,min}$ (valence band minimum energy), 16 atoms and 399 bands. The upper and lower plots show the cost to compute the valence and the conduction band contributions to ${\mathrm{\Sigma }}^{\left(\mathrm{model}\right)}\left(\omega \right)$, respectively. The position of the two minima are $(N_{v_w}=2,N_{p_w}=7)$, upper, and $(N_{c_w}=1,N_{p_w}=3)$, lower.