Low-scaling \textit{GW} calculation of quasi-particle energies within numerical atomic orbital framework

This paper presents a low-scaling space-time $GW$ algorithm within the numerical atomic orbital framework that leverages the localized resolution of identity technique to reduce computational complexity to $O(N^2)$ or better, enabling efficient and accurate quasi-particle energy calculations for systems with fewer than 100 atoms.

The Gist

The Big Picture: Predicting the Future of Materials

Imagine you are an architect trying to design a new solar panel or a super-fast computer chip. To do this, you need to know exactly how electrons (the tiny particles that carry electricity) behave inside the material.

Scientists use a powerful mathematical tool called GW to predict this behavior. It's like having a crystal ball that tells you exactly how much energy an electron needs to jump from a "sitting" state to a "running" state. This is crucial for understanding why silicon is a semiconductor and why copper is a conductor.

The Problem: The crystal ball is heavy.
The traditional way to use this tool (called the "canonical method") is incredibly slow. It's like trying to count every single grain of sand on a beach by picking them up one by one. If you want to study a small rock, it takes a few hours. If you want to study a mountain (a large molecule or a complex material), it would take your computer centuries to finish. This is because the time it takes grows exponentially: double the size of the system, and the time needed quadruples (or even more).

The Solution: A New "Space-Time" Shortcut

The authors of this paper (Zhang, Lin, Shi, and Ren) have built a new, much faster version of this crystal ball. They call it a "Low-Scaling GW calculation."

Here is how they did it, using three simple analogies:

1. The "Neighborhood" Strategy (Local Resolution of Identity)

The Old Way: Imagine you are trying to calculate the noise level in a massive city. The old method assumes that every person in the city is talking to every other person in the city simultaneously. You have to check 1 billion conversations. That's impossible.

The New Way: The authors realized that in the real world, people mostly talk to their neighbors. You don't need to check if someone in Tokyo is talking to someone in New York; they are too far apart to hear each other.
They use a technique called LRI (Localized Resolution of Identity). It's like saying, "We only need to calculate the conversations between neighbors." By ignoring the distant, silent connections, they cut the amount of work down from "checking the whole city" to "checking the neighborhood."

2. The "Real-Time" vs. "Frequency" Switch

The Old Way: Imagine trying to understand a song by looking at a sheet of music (the frequency domain). You have to analyze every note, every harmony, and every rhythm all at once. It's complex and requires a huge library of data.

The New Way: The authors use a Space-Time algorithm. Instead of looking at the whole song at once, they listen to the music second-by-second (the time domain).
They use a special trick (Fast Fourier Transforms) to switch back and forth between "listening to the moment" and "analyzing the notes." This allows them to process the data in a way that is much more efficient, turning a massive, tangled knot of math into a straight, easy-to-walk path.

3. The "Sparsity" Filter (Ignoring the Ghosts)

The Old Way: Even with the neighborhood strategy, the computer still tries to calculate tiny, invisible interactions that are so weak they don't matter. It's like trying to measure the weight of a feather while standing on a scale that also weighs the air around you.

The New Way: The authors added a filter. They set a rule: "If an interaction is weaker than a whisper, ignore it."
In the computer code, this is called matrix filtering. It's like a bouncer at a club who only lets in the important guests. By throwing out the "ghost" data that doesn't change the result, the computer runs lightning fast without losing accuracy.

The Results: From "Impossible" to "Doable"

The team tested their new method on real materials like Silicon (used in chips) and Diamond.

• Accuracy: They compared their new "neighborhood" method against the old "count-every-grain" method. The results were almost identical. The new method predicted the energy levels of electrons with the same precision as the old one, just much faster.
• Speed:
  • For small systems (like a tiny molecule), the old method was still okay.
  • But for larger systems (around 100 atoms), the new method became the clear winner.
  • For very large systems (hundreds of atoms), the new method was orders of magnitude faster. While the old method would take days or weeks, the new method could do it in hours.
• Scaling: The old method's time went up like a rocket ($N^4$). The new method's time went up like a gentle hill ($N^2$ or $N^{2.7}$). This means as you add more atoms, the new method barely breaks a sweat, while the old method gets overwhelmed.

Why This Matters

This paper is a breakthrough because it opens the door to studying real-world materials that were previously too big to simulate.

• Before: Scientists could only study perfect, tiny crystals.
• Now: They can study larger, more complex structures, defects in materials, and potentially new types of solar cells or batteries.

Think of it as upgrading from a bicycle to a high-speed train. You can now travel much further, much faster, to discover new materials that could power our future. The authors have essentially built a "fast lane" for quantum physics calculations.

Technical Summary

1. Problem Statement

The GW approximation within many-body perturbation theory is the gold standard for calculating accurate electronic band structures and quasi-particle energies in materials. However, its application to large-scale systems is severely limited by computational cost.

• Scaling Bottleneck: Conventional GW implementations (typically in reciprocal space with plane-wave bases) scale as $O(N^4)$ with system size ($N$), where the rate-limiting steps are the evaluation of the polarization function ($\chi_0$) and the self-energy ($\Sigma$).
• Existing Solutions: Space-time algorithms using plane waves have reduced scaling to $O(N^3)$, but they suffer from large prefactors and high memory requirements due to dense real-space grids.
• Gap: While Numerical Atomic Orbitals (NAOs) offer compact, localized basis sets that could exploit sparsity, a low-scaling space-time GW algorithm specifically optimized for the NAO framework had not yet been implemented.

2. Methodology

The authors developed a low-scaling space-time GW algorithm within the NAO framework, implemented in the LibRPA library and interfaced with the FHI-aims code. The core methodology relies on three pillars:

A. Space-Time Formalism in Imaginary Time

Instead of working in the frequency domain, the method operates in the imaginary time domain.

• The Green's function ($G$) and screened Coulomb interaction ($W$) are computed in real space and imaginary time.
• The self-energy is calculated as a simple product: $\Sigma(i\tau) = -G(i\tau)W(i\tau)$.
• Results are transformed to the frequency domain via Fourier transforms and analytically continued to the real axis using Padé approximants.

B. Localized Resolution of Identity (LRI)

To handle the four-center integrals inherent in GW calculations efficiently, the authors employ the LRI technique (also known as Pair-Atom Density Fitting).

• Mechanism: Products of NAO basis functions are expanded using a set of auxiliary basis functions (ABFs).
• Localization: Unlike global RI, LRI restricts the expansion of a product $\phi_i \phi_j$ to auxiliary functions centered only on the atoms $I$ and $J$ themselves (or their immediate neighbors).
• Impact: This drastically reduces the number of triple coefficients to compute and store, enabling the formulation of algorithms that scale better than $O(N^4)$.

C. Real-Space Tensor Contraction & Sparsity

The core computational engine utilizes real-space tensor contractions over atom-pair blocks.

• Block Sparsity: The algorithm exploits the spatial locality of NAOs. Matrices and tensors are treated as blocks of atom pairs.
• Filtering: A key innovation is the application of threshold-based filtering ($\eta$) on atom-pair blocks. Blocks with elements smaller than a threshold are discarded entirely before computation.
• Algorithmic Structure: The authors derived optimized three-layer loop structures for computing the response function ($\chi_0$) and the self-energy ($\Sigma$). By reordering loops to maximize the reuse of intermediate arrays and leveraging the sparsity of the LRI coefficients, they avoid unnecessary $O(N^4)$ operations.

D. Handling Singularities

• $\Gamma$-Point Singularity: The method addresses the divergence of the Coulomb matrix at the $\Gamma$ point ($q \to 0$) using a regularization scheme involving eigenvalue decomposition of the Coulomb matrix and applying corrections to the dielectric matrix (head and wing terms).
• Coulomb Truncation: An attenuated Coulomb kernel is used in real-space integrations to improve numerical stability.

3. Key Contributions

1. First Low-Scaling NAO-GW: Implementation of the first space-time GW algorithm specifically designed for the NAO basis set, achieving formal scaling better than $O(N^4)$.
2. Unified Tensor Framework: Development of a unified real-space tensor-contraction framework that handles both the polarization function and the self-energy efficiently.
3. Block-Wise Filtering: Introduction of a robust atom-pair block filtering strategy that significantly accelerates calculations without sacrificing accuracy, particularly for the correlation self-energy.
4. LibRPA Integration: Successful integration into the open-source LibRPA library, making these capabilities accessible to the FHI-aims ecosystem.

4. Results and Benchmarks

The authors validated the implementation against the canonical $O(N^4)$ $k$-space implementation in FHI-aims using a wide range of semiconductors and insulators (e.g., Si, MgO, Diamond, GaAs).

• Accuracy:

  • Band Gaps: The low-scaling implementation reproduces canonical results with high fidelity. With 32 minimax grid points, the mean absolute error (MAE) in fundamental band gaps is < 5 meV for most systems.
  • Band Structures: Quasi-particle band structures for Si and MgO are nearly indistinguishable from the canonical results up to 30 eV above the valence band maximum.

• Scaling with System Size ($N$):

  • $\chi_0$ (Response Function): Scales as $O(N^{2.5})$ (empirically), crossing over to be faster than the canonical method at ~80 atoms.
  • $\Sigma_c$ (Correlation Self-Energy): Scales as $O(N^{2.3})$, crossing over at ~40 atoms.
  • Overall Scaling: The total wall time scales as $O(N^{2.7})$. The $O(N^3)$ bottleneck shifts to the construction of the screened Coulomb interaction ($W^c$), which involves dense matrix operations.
  • Threshold: The low-scaling method becomes advantageous for systems containing fewer than 100 atoms.

• Scaling with k-points ($N_k$):

  • Canonical: $O(N_k^2)$.
  • Low-scaling: Linear scaling ($O(N_k)$) for dense meshes ($N_k > 216$). The crossover point is between $7\times7\times7$ and $8\times8\times8$ k-grids.

• Parallel Efficiency:

  • The code demonstrates strong scaling up to ~10,000 CPU cores (12,288 cores tested) for systems with 512 atoms, with parallel efficiencies around 20%.

5. Significance

This work represents a major step forward in making high-accuracy GW calculations feasible for large-scale materials (hundreds of atoms).

• Accessibility: It lowers the computational barrier, allowing researchers to study complex defects, interfaces, and nanostructures that were previously inaccessible to GW due to $O(N^4)$ costs.
• Efficiency: By combining the compactness of NAOs with the efficiency of the space-time method and LRI, the authors have achieved a scaling regime ($O(N^{2.7})$) that is significantly superior to traditional approaches.
• Future Potential: The framework is general and can be extended to other localized basis sets. Future work will focus on optimizing the $W^c$ construction (potentially via GPU acceleration), handling spin-orbit coupling, and extending the method to metallic systems and low-dimensional materials.