Implementation of the hybrid exchange-correlation functionals in the SIESTA code

This paper presents an efficient and accurate implementation of hybrid exchange-correlation functionals in the SIESTA code, utilizing a Gaussian-fitted representation of numerical atomic orbitals to enable large-scale, scalable simulations of extended systems with significantly improved band gap predictions.

The Gist

Imagine you are trying to simulate how a city of atoms behaves. You want to know how electrons (the tiny particles that hold atoms together) move and interact. For decades, scientists have used a tool called Density Functional Theory (DFT) to do this. Think of DFT as a very fast, very efficient map. It's great for getting a general idea of the city's layout, but it has a blind spot: it often gets the "energy gaps" (the distance between the ground floor and the first floor of a building) wrong. It tends to say the gap is smaller than it really is, which can make a material look like a conductor when it's actually an insulator.

To fix this, scientists developed Hybrid Functionals. These are like upgrading your map to include a high-definition satellite view. They add a specific type of "exact exchange" calculation that corrects the blind spots, giving you the right energy gaps. However, there's a catch: this high-definition view is incredibly slow to compute. It's like trying to calculate the traffic flow for every single car in a massive city simultaneously; the computer gets overwhelmed and the simulation takes forever.

The Problem: The "Four-Center" Bottleneck
The main reason hybrid calculations are so slow is a math problem involving "four-center integrals." Imagine trying to calculate the interaction between four different people in a room. If you have 1,000 people, the number of possible four-person groups is astronomical. In the world of atoms, calculating these interactions for every possible group is the computational bottleneck.

The Solution: The "Gaussian" Translator
The authors of this paper, working with the SIESTA code (a popular software for simulating materials), found a clever way to speed this up.

1. The Native Language (NAOs): SIESTA usually speaks in "Numerical Atomic Orbitals" (NAOs). These are like strict, localized maps that stop abruptly at a certain distance. They are efficient for standard calculations but very hard to use for the complex "four-center" math required for hybrid functionals.
2. The Translation (GTOs): The team created a translator. They took those strict, localized maps (NAOs) and approximated them using "Gaussian-type orbitals" (GTOs). Think of GTOs as smooth, bell-curve shapes that are mathematically friendly.
3. The Library (Libint): Because GTOs are mathematically smooth, there is a pre-existing, highly optimized "library" (called libint) that can instantly calculate the interactions between them. It's like having a pre-calculated dictionary for every possible conversation between four people.

How They Made It Work
The team didn't just swap the languages; they built a bridge:

• Fitting: They mathematically "fitted" the strict SIESTA maps into the smooth Gaussian shapes. It's like taking a jagged, pixelated image and smoothing it out so a high-end printer can handle it, without losing the original picture's details.
• Screening: They added a "bouncer" at the door. Since most atoms are too far apart to interact significantly, the code ignores those distant pairs. This reduces the number of calculations from billions to a manageable few million.
• Parallel Power: They built a system where thousands of computer processors can work on different parts of the city simultaneously without stepping on each other's toes.

The Results: Faster and More Accurate
The paper tested this new method on a wide variety of materials, from silicon chips to 2D materials like graphene.

• Accuracy: The new method fixed the "blind spots." For example, it correctly predicted that black phosphorus is a semiconductor (with a gap) rather than a metal, and it calculated the energy gaps of silicon and diamond to be almost identical to experimental reality.
• Speed: By using the Gaussian translation and the screening "bouncer," they made these high-accuracy calculations feasible for large systems (hundreds or even thousands of atoms) that would have previously taken too long to run.

The Trade-Off
The authors also analyzed how to get the best balance between speed and accuracy. They found that:

• Using a moderate number of "Gaussian shapes" (about 4 to 6) to represent each atom is usually enough.
• Setting a specific "cutoff" distance for interactions works well without needing to calculate every single distant atom.
• This balance allows scientists to get results that are nearly as accurate as the most expensive methods, but in a fraction of the time.

In Summary
This paper presents a new engine for the SIESTA software. It allows scientists to run high-precision, "hybrid" simulations on large materials by translating the software's native language into a mathematically smoother one that can be processed instantly. This makes it possible to accurately predict the electronic properties of complex materials (like semiconductors and 2D sheets) without waiting weeks for the computer to finish the job.

Technical Summary

1. Problem Statement

Density Functional Theory (DFT) is a cornerstone of materials science, but standard semi-local functionals (like LDA and GGA/PBE) suffer from the self-interaction error. This leads to systematic underestimation of band gaps, overestimation of electron delocalization, and inaccuracies in describing defect levels and reaction barriers.
Hybrid functionals (e.g., HSE06, PBE0), which mix a fraction of exact Hartree-Fock (HF) exchange with semi-local exchange, solve these issues. However, their application to large-scale extended systems (solids, surfaces, defects) within the SIESTA code has been historically limited due to:

• Computational Cost: The evaluation of four-center electron repulsion integrals (ERIs) required for exact exchange scales poorly ($O(N^4)$) and is computationally prohibitive.
• Basis Set Mismatch: SIESTA uses strictly localized Numerical Atomic Orbitals (NAOs), which offer excellent spatial locality and sparse matrix structures. However, standard efficient libraries for analytical ERI evaluation (like libint) are designed for Gaussian-Type Orbitals (GTOs).
• Lack of Analytical Forces: Previous implementations often lacked fully analytical force calculations, hindering geometry optimization and molecular dynamics.

2. Methodology

The authors present a hybrid scheme that leverages the strengths of both NAOs and GTOs to enable efficient hybrid functional calculations in SIESTA.

A. NAO-to-GTO Fitting Strategy

Instead of evaluating ERIs numerically (which is slow) or abandoning the NAO basis, the authors approximate the radial part of the native NAOs as a linear combination of Gaussians:
$$\chi_\mu(r) \approx \sum_{k=1}^{N_G} C_k \exp(-\alpha_k r^2)$$

• Internal Fitting: A new, robust fitting algorithm is implemented directly within SIESTA. It uses an alternating optimization strategy (fixing exponents to solve for coefficients via least squares, then optimizing exponents) to minimize the difference between the NAO and the Gaussian expansion.
• Strict Localization: Since Gaussians are non-zero at infinity, the authors introduce a truncation strategy. The Gaussian sum is multiplied by a smooth cutoff function ($1 - e^{-\alpha(r-r_c)^2}$) to ensure the orbitals remain strictly localized and continuous at the cutoff radius, preserving SIESTA's sparse matrix advantages.
• Angular Transformation: The real spherical harmonics used in SIESTA are transformed into Cartesian Gaussian harmonics to interface with the libint library.

B. Integral Evaluation and Screening

• Analytical Evaluation: Once NAOs are represented as contracted Gaussians, the four-center ERIs are computed analytically using the libint library (version 1.1.4), which employs the Obara-Saika algorithm.
• Multi-Level Screening: To reduce the computational complexity from $O(N^4)$ to near-linear scaling, the authors employ a hierarchy of screening techniques:
  1. Schwarz Inequality: Bounds the magnitude of integrals to discard negligible pairs.
  2. Shell Quartet Filtering: Groups orbitals into shells and filters based on spatial overlap and distance.
  3. Density Matrix Screening: Multiplies the Schwarz bound by the maximum density matrix element connecting the shells; if the product is below a threshold, the integral is skipped.
  4. Range-Separation: Utilizes the short-range nature of the HSE06 functional (screened Coulomb potential) to naturally truncate long-range interactions, eliminating the need for Ewald summation.

C. Analytical Forces and Parallelization

• Forces: The implementation derives fully analytical forces, including Pulay terms arising from the basis set dependence on atomic positions. This allows for geometry optimization and ab initio molecular dynamics.
• Parallelization: A dynamic parallel distribution scheme is used to distribute the calculation of ERIs across MPI processes, ensuring excellent load balancing and scalability.

3. Key Contributions

1. First Efficient Hybrid Implementation in SIESTA: Successfully integrates exact exchange into the NAO-based SIESTA framework without sacrificing the code's inherent linear-scaling capabilities.
2. Robust Internal Fitting Algorithm: Developed a self-contained, systematic method to fit NAOs to Gaussians, removing the need for external fitting tools and complex workflows.
3. Analytical Forces: Provided the first fully analytical force formulation for hybrid functionals within the SIESTA NAO framework, enabling structural relaxations and dynamics.
4. Comprehensive Benchmarking: Validated the implementation across a wide range of materials (semiconductors, insulators, 2D materials, and molecules) against PBE, $G_0W_0$, and experimental data.

4. Results

The paper validates the implementation through extensive benchmarks:

• Band Gaps: The HSE06 functional significantly improves band gap predictions compared to PBE.
  • Example: Silicon band gap improved from 0.59 eV (PBE) to 1.20 eV (HSE06), matching the experimental value of 1.12 eV.
  • Example: Black Phosphorus, which PBE incorrectly predicts as metallic, is correctly predicted as a semiconductor with a gap of 0.29 eV (close to $G_0W_0$ and experiment).
  • The results show excellent agreement with $G_0W_0$ calculations and experimental data across Si, Diamond, LiF, c-BN, TiO$_2$, ZnO, and h-BN.
• Accuracy vs. Efficiency Trade-offs:
  • Number of Gaussians ($N_G$): Increasing $N_G$ improves accuracy but drastically increases cost ($O(N_G^4)$). The authors find that $N_G = 4$ offers an optimal balance, providing gaps within tens of meV of the fully converged value with a fraction of the computational time.
  • Cutoff Radius ($\epsilon_{Gauss}$): A threshold of $5 \times 10^{-3}$ provides sufficiently compact orbitals to maintain sparsity while ensuring converged band gaps for most systems.
  • Screening Thresholds: Thresholds of $10^{-5}$ to $10^{-6}$ for integral screening provide nearly converged results with a significant speedup (order of magnitude) compared to unfiltered calculations.
• Parallel Scalability: Benchmarks on a 5$\times$5$\times$5 SrTiO$_3$ supercell (hundreds of atoms) show an 11.8$\times$ speedup when scaling from 4 to 64 cores, with parallel efficiency remaining above 73%.
• Molecular Validation: The method correctly reproduces the singlet-triplet splitting in the O$_2$ molecule, confirming the accuracy of the spin-resolved exchange implementation.

5. Significance

This work removes a major bottleneck in first-principles simulations of large extended systems. By enabling hybrid functional calculations for systems containing several hundred to a few thousand atoms within SIESTA, the authors have made accurate predictions of electronic structures (band gaps, defect levels) and structural properties accessible at a scale previously reserved for semi-local functionals or much smaller systems.

The ability to perform geometry optimizations and molecular dynamics with hybrid functionals opens new avenues for studying:

• Polaronic behavior and charge localization.
• Defect physics in semiconductors and insulators.
• Surface and interface properties in heterostructures.
• Accurate thermochemical properties for complex materials.

The implementation positions SIESTA as a competitive, versatile tool for high-accuracy materials modeling, bridging the gap between the efficiency of localized basis sets and the accuracy of hybrid functionals.