Local Basis DFT in SIESTA

• Local Basis DFT is a method that employs strictly confined numerical atomic orbitals for efficient and accurate electronic structure simulations.
• SIESTA optimizes basis functions through parameter tuning and soft confinement, achieving near-quantitative agreement with plane-wave benchmarks.
• The approach integrates advanced solvers like PEXSI to deliver scalable performance, making it ideal for low-dimensional materials like graphene and hBN.

Local-basis density functional theory (DFT) as realized in the SIESTA code represents a paradigm in first-principles electronic structure simulation, defined by the use of strictly confined numerical atomic orbitals (NAOs) as a basis set. By enforcing locality, SIESTA achieves computational scaling and efficiency advantages over plane-wave (PW) approaches, especially for large or low-dimensional systems. The construction, optimization, and validation of these NAO bases are essential to maintain accuracy across a variety of environments, motivating algorithmic advances and systematic benchmarking against PW methodologies (Bennett et al., 2024).

1. Numerical Atomic Orbitals and Basis Construction

In SIESTA, each basis function $\varphi_{\alpha nlm}(\mathbf{r})$ is written in a separable form as $$\varphi_{\alpha nlm}(\mathbf{r}) = R_{\alpha nl}(r) Y_{lm}(\Omega)$$, combining a numerically tabulated radial function $R_{\alpha nl}(r)$ and real spherical harmonics $Y_{lm}(\Omega)$. Strict locality is enforced by truncating $R_{\alpha nl}(r)$ to zero beyond a species- and channel-specific cutoff radius $r_c$. Smoothness of the orbitals is guaranteed via a soft confinement potential,

$$V_{\rm conf}(r) = V_0\,\frac{\exp\left[-\frac{r_c - r_i}{r - r_i}\right]}{r_c - r},$$

where $r_i \lesssim r_c$ and $V_0 \gg 1$ Ry are parameters variationally tuned for each basis shell.

Zeta-levels (SZ, DZ, TZ, QZ), corresponding to the number of radial functions per angular momentum, augment the basis flexibility; for example, DZPF denotes a basis with double-zeta valence, polarization $d$ (l+1), and additional $f$ (l+2) shells. Polarization shells (up to $l+2$) are generated by solving the confined atomic Schrödinger equation for each angular momentum in the presence of the same pseudopotential (Bennett et al., 2024).

2. Basis-set Optimization Strategies

Basis functions' performance depends critically on their spatial range and completeness. To systematically construct optimized bases, basis parameters $\{r_c, V_0, r_i\}$ and a notional atomic charge $Q$ are adjusted to minimize a composite enthalpy,

$H = E_{\rm KS} + p_{\rm basis}\,V_{\rm orbs},$

where $E_{\rm KS}$ is the self-consistent Kohn–Sham energy at a fixed, PW-converged lattice constant, $V_{\rm orbs} \propto r_c^3$ is the total basis volume, and $p_{\rm basis}$ is a fictitious basis pressure controlling orbital extent (e.g., $0.03$ GPa in the optimized sets versus $0.2$ GPa in native SIESTA). The simplex algorithm iteratively refines these parameters per shell—keeping the system geometry fixed—for a compact yet accurate representation (Bennett et al., 2024).

Polarization flexibility is crucial for two-dimensional materials: inclusion of $l+2$ ($f$) shells substantially improves the basis' angular completeness for systems with hexagonal symmetry such as graphene and hBN.

3. Benchmarking and Accuracy

Key physical quantities directly used for benchmark validation include:

• Total energy difference:

$$\Delta E = E_{\rm SIESTA}(a_{\rm pw}) - E_{\rm PW}(a_{\rm pw})$$

• Lattice constant strain:

$$\eta = \frac{a_{\rm SIESTA} - a_{\rm pw}}{a_{\rm pw}}$$

• Cohesive energy:

$E_{\rm coh} = E_{\rm mono} - \sum_i E^i_{\rm atom}$

Applying these metrics to monolayer graphene (PW reference: $a = 2.466$ Å), native DZP basis yields $\Delta E \approx +0.8$ eV/atom and $\eta \approx +0.5\%$, while an optimized DZPF basis achieves $\Delta E < 0.02$ eV/atom and $\eta < 0.05\%$; cohesive energy errors similarly drop below $0.03$ eV relative to PW. hBN shows analogous trends. The addition of 4$f$ polarization is particularly effective in reducing energy and structural discrepancies (Bennett et al., 2024).

4. Computational Performance and Scaling

SIESTA's local-basis approach yields favorable computational scaling, especially for low-dimensional or large real-space supercells. For a 2-atom unit cell (20×20×1 $k$-grid), optimized DZP bases show $\sim2\times$ speedup over ABINIT at 1000 eV cutoff. For $\Gamma$-only supercells (50 to $\sim$300 atoms), wall-time for real-space integrals remains roughly linear in system size ($O(N)$), whereas ABINIT's cost exceeds $O(N^2)$ as FFT grid sizes grow. On a 12×12 supercell, SIESTA can be over $10\times$ faster than ABINIT. Furthermore, optimized bases yield additional $10–20\%$ speedups due to reduced orbital extent compared to native SIESTA defaults (Bennett et al., 2024).

5. Integration with Pseudopotentials and XC Functionals

SIESTA employs norm-conserving ONCVPSP pseudopotentials (Pseudo-Dojo library), read via the PSML interface, which ensures consistent local and nonlocal (Kleinman–Bylander) projectors between codes. The exchange-correlation potential is typically treated at the PBE GGA level through the LibXC library, and a real-space mesh cutoff of $1000$ Ry is used to represent densities and potentials—matching or exceeding the convergence standards of plane-wave DFT calculations (e.g., $>2000$ eV cutoff in ABINIT) (Bennett et al., 2024).

6. Advanced Solvers and O(N) Methods

To exploit matrix sparsity inherent to localized bases, SIESTA supports solvers providing sub-cubic scaling:

• PEXSI (Pole Expansion and Selected Inversion): Achieves $O(N^{1.5})$ scaling for 2D systems, $O(N^2)$ for 3D, and $O(N)$ for 1D (Lin et al., 2014). PEXSI does not require eigenvalue computation and computes physical observables through pole expansions of the Fermi operator and selective inversion for needed matrix elements.
• Nested Dissection and O($N^{7/3}$) methods: These achieve lower-order scaling for density matrix construction using contour integration of the Green function with recursive blockwise sparse inversion (Ozaki, 2010). Both approaches are numerically exact up to discretization/tolerance errors, with accuracy matching conventional cubic diagonalization, and are readily integrated into SIESTA's infrastructure.

Parallel scalability is leveraged by distributing either poles or domain hierarchies across distributed memory, enabling simulations of systems with dimensions $N$ up to $10^5$ orbitals on thousands of cores (Lin et al., 2014, Ozaki, 2010).

7. Application Scope and Systematic Validation

Local-basis DFT in SIESTA is widely used for low-dimensional materials (e.g., graphene, hBN, twisted bilayer graphene), interfaces, molecular adsorption, and superconductors. Optimized NAO bases confer near-quantitative agreement with plane-wave reference data for both atomic and electronic structures. This is evidenced in large moiré supercell studies of twisted bilayer graphene, where SIESTA with optimized bases reproduces continuum model trends and ab initio benchmarks, subject to systematic, quantifiable basis (and functional/pseudopotential) effects (Zhu et al., 23 Jan 2026). In molecular and surface science, explicit protocols for basis optimization and superposition error correction (BSSE) have been established, allowing direct comparison to experiment and PW-level accuracy (Buimaga-Iarinca et al., 2014). The local basis methodology is continuously being extended (TDDFT, superconducting DFT, transport, and perturbation theory) within SIESTA's framework (García et al., 2020, Reho et al., 2024, Coulaud et al., 2013).

---

In conclusion, the local-basis DFT formalism as implemented in SIESTA, when equipped with carefully optimized and polarized NAO sets, achieves rapid, scalable, and accurate electronic structure simulations rivaling plane-wave benchmarks, particularly excelling in large-scale and low-dimensional material applications. Incorporation of advanced sparse solvers and ensured interoperability with modern pseudopotential and XC frameworks cements its role in first-principles materials modeling (Bennett et al., 2024, Ozaki, 2010, Lin et al., 2014, García et al., 2020, Zhu et al., 23 Jan 2026, Buimaga-Iarinca et al., 2014, Coulaud et al., 2013, Reho et al., 2024).