The SIESTA method for ab initio order-N materials simulation

TL;DR

The paper presents Siesta, an ab initio DFT method designed for order-N scaling in large materials systems by combining a flexible localized LCAO basis with a real-space grid, and by recasting the Kohn-Sham problem using a localization-based energy functional and Wannier-like states. It provides full treatment of norm-conserving pseudopotentials (KB form with NLCC), multiple-zeta basis sets (SZ/DZP), and grid-based evaluation of Hartree and XC potentials, enabling efficient two-center and grid integrals and linear-scaling density updates. The implementation supports non-collinear spin, Brillouin-zone sampling via an auxiliary supercell approach, and optional finite-temperature occupations, while offering extensive features for structural relaxation, MD, phonon analysis, transport, and TDDFT. Collectively, Siesta delivers accurate, scalable first-principles simulations for thousands of atoms, with demonstrated convergence and computational efficiency relative to plane-wave approaches.

Abstract

We have developed and implemented a self-consistent density functional method using standard norm-conserving pseudopotentials and a flexible, numerical LCAO basis set, which includes multiple-zeta and polarization orbitals. Exchange and correlation are treated with the local spin density or generalized gradient approximations. The basis functions and the electron density are projected on a real-space grid, in order to calculate the Hartree and exchange-correlation potentials and matrix elements, with a number of operations that scales linearly with the size of the system. We use a modified energy functional, whose minimization produces orthogonal wavefunctions and the same energy and density as the Kohn-Sham energy functional, without the need of an explicit orthogonalization. Additionally, using localized Wannier-like electron wavefunctions allows the computation time and memory, required to minimize the energy, to also scale linearly with the size of the system. Forces and stresses are also calculated efficiently and accurately, thus allowing structural relaxation and molecular dynamics simulations.

Figures (8)

• Figure 1: Local pseudopotential for silicon. $V_{local}$ is the unscreened local part of the pseudopotential, generated as the electrostatic potential produced by a localized distribution of positive charge, Eq. (\ref{['rho_local']}), whose integral is equal to the valence ion charge ($Z=4$ for Si). The dashed line is $-Z/r$. $V_{NA}$ is the local pseudopotential screened by an electron charge distribution, generated by filling the first-$\zeta$ basis orbitals with the free-atom valence occupations. Since these basis orbitals are strictly confined to a radius $r^c_{max}$, $V_{NA}$ is also strictly zero beyond that radius.
• Figure 2: Comparison of convergence of the total energy with respect to the sizes of a plane wave basis set and of the LCAO basis set used by Siesta. The curve shows the total energy per atom of silicon versus the cutoff of a plane wave basis, calculated with a program independent of Siesta, which uses the same pseudopotential. The arrows indicate the energies obtained with different LCAO basis sets, calculated with Siesta, and the plane wave cutoffs that yield the same energies. The numbers in parentheses indicate the basis sizes, i.e. the number of atomic orbitals or plane waves of each basis set. SZ: single-$\zeta$ (valence $s$ and $p$ orbitals); DZ: double-$\zeta$; TZ: triple-$\zeta$; DZP: double-$\zeta$ valence orbitals plus single-$\zeta$ polarization $d$ orbitals; TZP: triple-$\zeta$ valence plus single-$\zeta$ polarization; TZDP: triple-$\zeta$ valence plus double-$\zeta$ polarization; TZTP: triple-$\zeta$ valence plus triple-$\zeta$ polarization; TZTPF: same as TZTP plus extra single-$\zeta$ polarization $f$ orbitals.
• Figure 3: Total energy per atom versus lattice constant for bulk silicon, using different basis sets, noted as in Fig. \ref{['PWcomparison']}. PW refers to a very well converged (50 Ry cutoff) plane wave calculation. The dotted line joins the minima of the different curves.
• Figure 4: Dependence of the lattice constant, bulk modulus, and cohesive energy of bulk silicon with the cutoff radius of the basis orbitals. The $s$ and $p$ orbital radii have been made equal in this case, to simplify the plot. PW refers to a well converged plane wave calculation with the same pseudopotential.
• Figure 5: (a) Convergence of the total energy and pressure in bulk silicon as a function of the energy cutoff $E_{cut}$ of the real space integration mesh. Circles and continuous line: using a grid-cell-sampling of eight refinement points per original grid point. The refinement points are used only in the final calculation, not during the self-consistency iteration (see text). Triangles: two refinement points per original grid point. White circles: no grid-cell-sampling. (b) Bond length and angle of the water molecule as a function of $E_{cut}$