Efficient all-electron Bethe-Salpeter implementation using crystal symmetries

This paper presents an efficient all-electron Bethe-Salpeter equation implementation within the FLAPW method that leverages crystal symmetries to block-diagonalize the two-particle Hamiltonian, significantly accelerating diagonalization and yielding improved exciton binding energies for materials like Si, LiF, and MoS$_2$.

The Gist

Imagine you are trying to predict how a crystal, like a piece of silicon or a flake of molybdenum disulfide, will react when light hits it. You want to know exactly which colors it absorbs and which it reflects. To do this, scientists use a complex mathematical recipe called the Bethe-Salpeter Equation (BSE).

Think of the BSE as a giant, incredibly detailed simulation of a dance floor inside the crystal. On this floor, electrons (the dancers) and "holes" (the empty spaces they leave behind) pair up to form "excitons" (dance couples). The goal is to calculate the energy of these couples to predict the crystal's optical properties.

However, there are two major problems with this dance floor simulation:

1. The Crowd is Too Big: To get an accurate result, you need to simulate millions of dancers on a grid that covers the entire crystal. The math required to track every single pair becomes so massive that even the world's fastest supercomputers struggle to solve it.
2. The "All-Electron" Challenge: Most previous simulations took a shortcut. They pretended the inner core of the atoms didn't exist, focusing only on the outer dancers. This is like trying to choreograph a dance while ignoring the dancers' legs. The authors of this paper wanted to do the full, "all-electron" version, which is much more accurate but computationally terrifying.

The Solution: The Crystal's Secret Symmetry

The authors, working at RWTH Aachen and Jülich, found a brilliant way to speed this up without losing accuracy. They realized that crystals are not random messes; they are built with perfect symmetry. If you rotate a crystal or flip it, it looks exactly the same.

They used this symmetry like a magic key to unlock the problem in three clever ways:

1. The "Copy-Paste" Trick (Building the Hamiltonian)

Imagine you are filling out a massive spreadsheet where every cell represents a possible interaction between two dancers. Usually, you'd have to calculate the value for every single cell.

• The Old Way: Calculate 1,000,000 cells one by one.
• The New Way: Because the crystal is symmetrical, if you calculate the interaction for one specific pair of dancers, you can mathematically deduce the interactions for thousands of other pairs just by rotating or flipping that result.
• The Result: They only had to calculate a tiny fraction of the spreadsheet explicitly. The rest was filled in automatically using symmetry rules.

2. The "VIP Lounge" (Diagonalization)

Once the spreadsheet (called the Hamiltonian matrix) is built, the computer has to solve it to find the energy levels. This is like trying to find a specific needle in a haystack the size of a mountain.

• The Old Way: Tackle the whole mountain of data at once.
• The New Way: The authors used group theory (a branch of math about symmetry) to rearrange the spreadsheet. They discovered that the massive mountain of data could be split into separate, smaller rooms (blocks).
• The Magic: They found that for light absorption, only one of these rooms actually matters. The other rooms are empty or irrelevant for this specific task.
• The Result: Instead of solving a problem with a million variables, they only had to solve a problem with a few hundred thousand. In the case of Silicon, this reduced the work by a factor of 125. It's like realizing you only need to check the VIP lounge for the party, not the entire stadium.

3. Handling the "Infinite" Problem

There was also a tricky mathematical issue where the calculations went to infinity (like dividing by zero) when electrons were very close together. The authors developed a sophisticated way to smooth out these infinities, treating them like a gentle curve rather than a sharp spike, ensuring the math stayed stable even for complex, layered materials like MoS2.

The Results: A Crystal Clear Picture

By using these symmetry tricks, the team successfully ran the "all-electron" simulation on three different materials:

• Silicon (Si): They achieved a level of detail never seen before, matching experimental data almost perfectly.
• Lithium Fluoride (LiF): A material with a huge energy gap, where the electron-hole pairs are very tightly bound.
• Molybdenum Disulfide (MoS2): A layered material used in advanced electronics.

Why does this matter?
Before this paper, scientists had to choose between accuracy (using the all-electron method) and speed (using the shortcut methods). This paper proved you can have both. By using the crystal's own symmetry as a shortcut, they made the impossible possible.

In a nutshell:
The authors built a super-accurate simulator for how crystals interact with light. Instead of brute-forcing the calculation, they realized the crystal's symmetry allowed them to ignore 80% of the work and focus only on the one part that actually matters. This turned a task that would take a supercomputer weeks into one that takes days, giving us a clearer, more accurate picture of the quantum world.

Technical Summary

1. Problem Statement

The Bethe-Salpeter Equation (BSE) is the standard first-principles method for calculating optical absorption spectra and exciton binding energies in solids. However, existing implementations face two major challenges:

• All-Electron Accuracy vs. Efficiency: Most high-accuracy methods (like FLAPW) resort to plane-wave basis sets for the interaction potentials (Coulomb and screened Coulomb) to handle singularities. This forces a "pseudopotential-like" approximation for the potentials near atomic nuclei, losing the true all-electron description of the wavefunctions in the core region.
• Computational Cost: Solving the BSE requires diagonalizing a large two-particle Hamiltonian matrix. The dimension of this matrix scales with the number of $k$-points ($N_k$), occupied bands ($N_o$), and unoccupied bands ($N_u$). Since accurate convergence requires dense $k$-point grids, the matrix becomes prohibitively large, making diagonalization the most time-consuming step (scaling cubically with matrix size).

2. Methodology

The authors present an implementation within the SPEX code (part of the FLEUR family) that addresses these issues through three main technical pillars:

A. All-Electron Mixed Basis for Interactions

Instead of projecting interaction potentials onto a plane-wave basis, the authors expand the bare ($v$) and screened ($W$) Coulomb potentials in a mixed basis.

• This basis is derived from the Linearized Augmented Plane-Wave (LAPW) method.
• In the interstitial region, it uses plane waves; inside the muffin-tin spheres, it uses products of radial functions and spherical harmonics.
• Significance: This preserves the full all-electron description of the wavefunctions and potentials, avoiding the smoothing artifacts associated with pseudopotentials or plane-wave projections near nuclei.
• Singularity Handling: The authors develop a rigorous method to handle the $1/q^2$ divergence of the screened interaction at $q \to 0$. They employ a modified integration technique (using an exponential cutoff and spherical harmonic expansion) that accounts for anisotropic dielectric screening (crucial for layered materials like MoS$_2$), avoiding the need for double-counting corrections that often plague other methods.

B. Exploitation of Crystal Symmetries

The core innovation is the systematic use of spatial and time-reversal symmetries to accelerate the calculation in two stages:

1. Hamiltonian Construction: Instead of calculating all matrix elements explicitly, the code calculates a minimal set of "seed" elements within the irreducible Brillouin zone (IBZ) and generates the rest via symmetry transformations.
2. Block-Diagonalization: Using group theory, the authors transform the electron-hole product basis into a symmetry-adapted basis.
   • This transforms the dense Hamiltonian matrix into a block-diagonal form, where each block corresponds to an irreducible representation (irrep) of the crystal's symmetry group.
   • Key Insight: For optical absorption, the oscillator strength is non-zero only for specific irreps. In many cases (including the examples studied), only one block contributes to the optical spectrum. This allows the diagonalization to be performed on a single, much smaller block rather than the full matrix.

C. Parallelization and Storage

To handle the massive memory requirements of large systems:

• The Hamiltonian is stored in distributed memory across multiple nodes using MPI.
• A specific strategy is employed to handle non-contiguous memory writes caused by symmetry operations.
• The transformation to the symmetry-adapted basis is performed in a way that produces a block-cyclic distribution, making the matrix immediately compatible with high-performance parallel diagonalization libraries (ScaLAPACK or ELPA).

3. Key Contributions

• First All-Electron FLAPW BSE with Mixed Basis: Demonstrates a fully all-electron implementation where interaction potentials are expanded in the mixed basis, retaining accuracy near atomic nuclei without plane-wave approximations.
• Group-Theoretical Acceleration: Introduces a method to reduce the BSE Hamiltonian to block-diagonal form. The paper shows that for optical spectra, often only one block is relevant, leading to massive speedups.
• Anisotropic Screening Treatment: Provides a robust mathematical framework for treating the Coulomb divergence in systems with non-diagonal dielectric tensors (anisotropic screening).
• Scalable Parallel Implementation: Develops efficient algorithms for distributing non-contiguous matrix elements and transforming them into block-cyclic formats suitable for exascale computing.

4. Results and Validation

The implementation was validated on three materials: Silicon (Si), Lithium Fluoride (LiF), and Molybdenum Disulfide (MoS$_2$).

• Silicon (Si):

  • Calculated on an extremely dense 60×60×60 $k$-point grid.
  • Symmetry reduced the Hamiltonian dimension by a factor of 5.
  • This resulted in a 125-fold speedup in the diagonalization step ($5^3$).
  • Exciton Binding Energy: Calculated as 22.5 meV, closer to the experimental value (~15 meV) than previous BSE studies.
  • The optical absorption spectrum matched experimental lineshapes very well.

• Lithium Fluoride (LiF):

  • Handled a Hamiltonian with nearly 2 million rows (full dimension).
  • Symmetry reduced the active block to 0.38 million rows.
  • Exciton Binding Energy: Extrapolated to 2.16 eV (optical), consistent with recent high-accuracy studies and experimental data.
  • Successfully reproduced the sharp 1s exciton peak at 12.6 eV.

• Molybdenum Disulfide (MoS$_2$):

  • Included Spin-Orbit Coupling (SOC) and treated anisotropic screening.
  • Symmetry reduced the dimension by a factor of 6, yielding a 216-fold speedup ($6^3$).
  • Exciton Binding Energy: Calculated as 76 meV, significantly closer to experimental values (~85 meV) than previous studies (which often reported ~130 meV due to coarse $k$-grids).
  • Correctly reproduced the characteristic double-peak structure (Excitons A and B) in the absorption spectrum.

5. Significance

This work represents a major advancement in computational materials science:

1. Accuracy: It enables truly all-electron BSE calculations, removing the approximations associated with pseudopotentials and plane-wave projections of potentials, which is critical for systems with strong core-valence interactions.
2. Efficiency: By leveraging crystal symmetries to reduce the diagonalization problem from $O(N^3)$ on a massive matrix to $O((N/r)^3)$ on a small block (where $r$ is the symmetry reduction factor), the method makes high-precision calculations on dense $k$-grids feasible.
3. Reliability: The results demonstrate that previous discrepancies in exciton binding energies (e.g., in MoS$_2$) were largely due to insufficient $k$-point sampling and methodological approximations, rather than fundamental flaws in the BSE formalism.
4. Scalability: The parallel implementation allows for the study of large, complex systems and the convergence of properties that were previously computationally inaccessible.

In summary, the paper provides a robust, efficient, and accurate framework for calculating optical properties of solids, bridging the gap between high-accuracy all-electron methods and the computational demands of many-body perturbation theory.