Source-linked AI summary
Daubechies wavelets as a basis set for density functional pseudopotential calculations
Luigi Genovese, Alexey Neelov, Stefan Goedecker, Thierry Deutsch, Seyed Alireza Ghasemi, Alexander Willand, Damien Caliste, Oded Zilberberg, Mark Rayson, Anders Bergman, Reinhold Schneider
TL;DR
The paper addresses the need for accurate, localized, and systematically convergent basis sets for isolated or inhomogeneous DFT calculations. It develops a Daubechies-wavelet pseudopotential method with analytic operators and efficient convolutions, implemented in ABINIT. The resulting code shows high systematic convergence, good performance, and parallel efficiency above 88%.
Problem
Plane-wave and non-systematic basis sets are poorly suited to localized or inhomogeneous systems and can introduce basis-choice uncertainty or numerical instability.
Method
The paper develops an isolated-system pseudopotential DFT method using a systematic orthogonal Daubechies-wavelet basis localized in real and Fourier space.
Results
Overall parallel efficiency is always higher than 88%, including for large systems with many processors.
Takeaways & Limitations
The method provides a freely available ABINIT implementation combining systematic convergence with efficient parallel electronic-structure calculations.
Abstract
from arXiv · showhide
Daubechies wavelets are a powerful systematic basis set for electronic structure calculations because they are orthogonal and localized both in real and Fourier space. We describe in detail how this basis set can be used to obtain a highly efficient and accurate method for density functional electronic structure calculations. An implementation of this method is available in the ABINIT free software package. This code shows high systematic convergence properties, very good performances and an excellent efficiency for parallel calculations.
I. INTRODUCTION
The introduction motivates a basis set that combines systematic convergence, orthogonality, real-space localization, Fourier-space localization, and adaptivity for isolated or inhomogeneous systems. It presents Daubechies wavelets as meeting these requirements and outlines their implementation in ABINIT.
- Motivation: Plane waves are inefficient for localized functions and isolated or inhomogeneous systems because they require large cells and high memory.Localized molecular charge densities cannot be exploited efficiently by a non-localized basis.
- Motivation: Systematic basis sets link numerical precision to the number of basis functions, reducing uncertainty from basis-set choice.Non-systematic Gaussian bases can become overcomplete and unstable before absolute convergence.
- Motivation: Nonorthogonal systematic bases require overlap-matrix operations, making them more complicated and slower.Orthogonality avoids these additional operations.
- Daubechies wavelets: Daubechies wavelets provide a systematic, orthogonal, smooth basis localized in both real and Fourier space while allowing adaptivity.Fourier localization supports preconditioning because high-frequency basis functions approximate kinetic-energy eigenfunctions.
- Method overview: Order-16 Daubechies wavelets are used; decreasing grid spacing increases the basis and improves numerical precision, with smoothness controlling convergence speed.The paper also describes linear-scaling fast wavelet transformations and more complex higher-level adaptivity.
- Daubechies wavelets: Compact support permits basis functions to be placed only near atoms, while pseudopotentials reduce the need for adaptivity to two resolution levels.The method uses high resolution around chemical bonds and lower resolution farther from atoms.
III. OVERVIEW OF THE METHOD
The method applies Kohn-Sham DFT to isolated systems using norm-conserving GTH-HGH pseudopotentials and a Daubechies scaling-function/wavelet representation. It supports analytic and tensor-product treatments of Hamiltonian terms, followed by direct minimization for gapped systems.
- Kohn-Sham formulation: The electronic density is computed from the squared moduli of Kohn-Sham wavefunctions, which are eigenfunctions of a pseudopotential Kohn-Sham Hamiltonian.The description assumes a closed-shell, nonspin-polarized system.
- Kohn-Sham formulation: The method targets isolated systems with free boundary conditions and includes Hartree, exchange-correlation, and external ionic potentials.The Hartree potential is obtained from Poisson’s equation.
- Pseudopotentials: It uses norm-conserving GTH-HGH pseudopotentials with local and nonlocal terms.Their analytic real-space form can be expressed using tensor products of one-dimensional functions.
- Optimization: After Hamiltonian application, Kohn-Sham wavefunctions are updated through direct minimization, implemented for non-zero-gap systems.The paper concentrates on insulators and notes possible extension to metals.
IV. TREATMENT OF KINETIC ENERGY
The kinetic-energy operator is represented analytically in the Daubechies basis and applied through short-range convolutions. This yields linear scaling with the number of nonzero wavefunction coefficients and an overall h14 convergence rate.
- Operator representation: Kinetic-energy matrix elements between scaling functions and wavelets are calculated analytically.The kinetic-energy expression is exact within a given Daubechies basis.
- Operator application: The projected kinetic-energy coefficients are related by a convolution using a three-dimensional filter formed from one-dimensional filters.This exploits the tensor-product structure of the basis.
- Operator application: Linear scaling is obtained for kinetic-energy evaluation with respect to the number of nonvanishing wavefunction expansion coefficients.For the order-16 family, the one-dimensional filter length is 29.
- Convergence: The wavefunction approximation error decreases as h8, while variational kinetic-energy error decreases as h14.The kinetic energy limits the overall convergence rate, which is shown in Figure 3.
V. TREATMENT OF LOCAL POTENTIAL ENERGY
The method avoids inaccurate individual matrix-element integration by smoothing Daubechies scaling functions and applying magic filters to obtain efficient local-potential expectation values. Its potential-energy evaluation converges as h^16, while the overall wavelet-code convergence is limited by the kinetic energy.
- Accurate local-potential evaluation is difficult in a Daubechies wavelet basis, motivating direct approximation of expectation values rather than individual matrix elements.The method uses matrix elements with respect to smoothed scaling functions.
- The magic filter maps scaling-function coefficients to smoothed functional values that accurately approximate local-potential expectation values.The approximation is not particularly accurate for a single matrix element but is excellent for expectation values.
- Wavelet and scaling-function operations are combined into modified magic filters on the original grid before applying the local potential.The transformation is needed because local-potential operations occur on a double-resolution grid with spacing h′ = h/2.
- The fine-scale filters combine a wavelet transform with the magic filter for Daubechies-16 scaling functions and wavelets.Figures 5 and 6 show the corresponding fine-scale filters, while Figure 4 shows the underlying magic filter.
- The local potential energy converges at h^16, two powers of h faster than the kinetic-energy convergence rate.
VI. CALCULATION OF HARTREE POTENTIAL
The charge density is constructed on a double-resolution grid and supplied to Poisson solvers using interpolating scaling functions. The resulting electrostatic treatment is accurate, efficient for large systems, and directly supports isolated charged systems.
- The discrete charge density on the double-resolution grid accurately represents the continuous charge distribution, with monopoles converging as h^16 and higher multipoles one power more slowly per order.Dipoles converge as h^15 and quadrupoles as h^14.
- For isolated molecules, interpolating scaling functions of order 16 conserve all multipoles through angular moment ℓ = 15 and impose correct free-boundary conditions.
- The Poisson solvers converge as h′^m; with order-16 interpolating scaling functions, electrostatic-potential convergence is faster than kinetic-energy convergence.
- Zero-padded FFT convolutions give the Poisson solvers an O(N log N) operation count.
- The Poisson solver requires roughly 1% of computational time for large systems and treats isolated net charges without compensating charges.Uniform potential accuracy allows use of the smallest volume compatible with decayed wavefunction tails at the boundary.
VII. XC FUNCTIONALS AND IMPLEMENTATION OF GGA’S
The implementation evaluates exchange-correlation terms from the real-space charge-density representation and supports both LDA and GGA functionals, including spin-polarised ABINIT routines. GGA gradients are computed on the double-resolution grid.
- Real-space charge-density data can be used directly to evaluate exchange-correlation energy and potential with ABINIT’s XC routines.The implementation also supports spin-polarised collinear ABINIT XC functionals.
- GGA functionals: GGA exchange-correlation energy density depends on both the local charge density ρ and the modulus of its gradient.
- GGA functionals: A fourth-order finite-difference scheme computes charge-density gradients on the double-resolution grid.
- GGA functionals: For free boundaries, charge-density values outside the computational volume are set equal to the corresponding border values when gradient stencils require them.
- GGA functionals: The White-Bird formulation accounts for the density-gradient dependence by splitting the exchange-correlation potential into ordinary and correction terms.The correction term appears only when the XC energy depends explicitly on |∇ρ|.
- XC evaluation and charge-density-gradient calculation can be combined with the Hartree Poisson solve to save computational time.
VIII. TREATMENT OF THE NON-LOCAL PSEUDOPOTENTIAL
Non-local pseudopotentials are evaluated efficiently by representing projectors in the same orthogonal Daubechies basis as the wavefunctions. The implementation uses separable Gaussian-polynomial projectors and direct minimisation with preconditioned iterative optimization.
- Non-local pseudopotential: Using the same Daubechies representation for projectors and wavefunctions reduces non-local pseudopotential operations to scalar products and coefficient updates.Orthogonality eliminates overlap-matrix operations in these projector applications.
- Non-local pseudopotential: GTH-HGH projectors, written as Gaussians multiplied by polynomials, are particularly convenient to expand in the Daubechies basis.
- Non-local pseudopotential: Separable three-dimensional projectors reduce the required integrals to products of three one-dimensional integrals.
- Non-local pseudopotential: Magic-filter quadrature on a grid 16 times denser gives an h^16 convergence rate for the one-dimensional projector integrations.The dense-grid procedure yields expansion coefficients accurate to machine precision for reasonable grid spacings.
- Minimisation and preconditioning: Direct total-energy minimisation uses preconditioned steepest descent or DIIS, with convergence defined by an average residual norm below a user-set tolerance.The gradient includes Lagrange multipliers enforcing wavefunction orthogonality.
- Minimisation and preconditioning: The preconditioned gradient is obtained by solving a discretized linear system with preconditioned conjugate gradients and a scaling-function-wavelet diagonal preconditioner.The initial solve uses typically four wavelet-resolution levels.
X. ORTHOGONALIZATION
The method maintains orthonormal Kohn–Sham wavefunctions during minimization by enforcing an identity overlap matrix. Orthogonalization dominates large-system costs, motivating a parallel-aware algorithm choice.
- X. ORTHOGONALIZATION: The wavefunctions must remain orthonormal throughout minimization, so the overlap matrix S must equal the identity matrix.This condition is the target of the orthogonalization step.
- X. ORTHOGONALIZATION: All orthogonalization algorithms have cubic complexity, making this stage dominant for large systems.The implementation therefore focuses on optimizing orthogonalization performance.
- A. Gram-Schmidt orthogonalization: On parallel computers, Gram–Schmidt requires n(n + 1)/2 communication steps when orbital coefficients are distributed, producing latency overhead and poor performance.The communication burden grows with the number of orbitals.
B. Loewdin orthogonalization
The paper compares Löwdin and Cholesky-based orthogonalization while also describing localized pseudopotential force evaluation. Cholesky avoids an eigenvalue problem and needs only one parallel communication step.
- B. Loewdin orthogonalization: Löwdin orthonormalization obtains new orbitals by multiplying the original set by the inverse square root of the overlap matrix S.This requires calculating S and solving an eigenvalue problem for its eigenvectors and eigenvalues.
- C. Pseudo Gram-Schmidt using Cholesky Factorization: The pseudo-Gram–Schmidt scheme factors the overlap matrix as S = LL^T before constructing the new orthonormal orbitals.The resulting orbitals are equivalent to those from classical Gram–Schmidt.
- C. Pseudo Gram-Schmidt using Cholesky Factorization: Cholesky factorization is faster than solving the Löwdin eigenvalue problem, giving a lower prefactor and requiring only one communication step in parallel.The overlap-matrix contribution assembly remains parallelized as in the Löwdin approach.
- Atomic forces: Atomic forces are evaluated directly through the Feynman–Hellmann theorem because the fixed basis produces no Pulay forces.Only the trivial ion–ion interaction is treated separately.
- Atomic forces: Localized pseudopotential force integrals can be restricted to small atom-centered regions and distributed across processors, with nearly linear O(N log N) scaling.This applies to the local pseudopotential contribution.
XII. LOCALIZATION PROPERTIES AND SMOOTHNESS OF THE BASIS FUNCTIONS
Localization parameters place basis functions only near atoms, using high resolution for bonds and low resolution farther away. Finite computational boundaries reduce tail smoothness and can affect kinetic-energy accuracy.
- XII. LOCALIZATION PROPERTIES AND SMOOTHNESS OF THE BASIS FUNCTIONS: Basis functions are associated with grid points inside unions of atom-centered spheres for both high- and low-resolution grids.The localization radii are selected to reduce degrees of freedom while meeting a target accuracy.
- XII. LOCALIZATION PROPERTIES AND SMOOTHNESS OF THE BASIS FUNCTIONS: The high-resolution region contains chemical bonds, while low-resolution regions represent exponentially decaying wavefunction tails.The localization parameters determine the sizes of these regions.
- XII. LOCALIZATION PROPERTIES AND SMOOTHNESS OF THE BASIS FUNCTIONS: Near a finite computational boundary, fewer scaling functions contribute, so wavefunctions become less smooth toward the edge.The reduced smoothness primarily affects the kinetic energy.
- XII. LOCALIZATION PROPERTIES AND SMOOTHNESS OF THE BASIS FUNCTIONS: For neutral systems, potential-energy errors decrease exponentially as the computational volume grows because the potential is very small far from atoms.This behavior concerns the finite-volume tail region.
XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS
Finite-volume truncation can create kinetic-energy errors in wavefunction tails, so the paper corrects these tails perturbatively after a medium-box self-consistent calculation. The correction targets a thin boundary shell at lower cost than a very large self-consistent volume.
- XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS: Tail kinetic-energy errors grow as h decreases when the computational volume is too small, because boundary nonsmoothness contributes approximately A/h^2 instead of A.The tail amplitude A is controlled by the wavefunction’s asymptotic decay and KS eigenvalue.
- XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS: With a sufficiently large computational volume, the method exhibits strict variational convergence at a rate of h^14 across a broad range of grid spacings.The convergence behavior is illustrated in Figure 3.
- XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS: The tail-correction method performs a self-consistent calculation in a medium box, then adds the missing far tail and cancels surface nonsmoothness afterward.The corrected wavefunction is represented as the medium-box wavefunction plus a perturbative correction.
- XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS: The finite-volume gradient is nonzero only in a small shell just outside the original volume, whose width is set by the kinetic-energy filter length.Farther outside the surface region, the gradient vanishes because the truncated wavefunction is identically zero.
- XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS: Figure 8 compares total-energy convergence versus low-resolution localization radius with and without tail corrections for two grid spacings.The plot tests h convergence once the localization parameter is sufficiently extended.
- XIII. PERTURBATIVE CALCULATION OF THE FINITE SIZE CORRECTIONS: The parallel implementation distributes orbitals or coefficients across processors, with Hamiltonian application and preconditioning mainly using orbital distribution.These schemes support parallel execution of the correction workflow.
XV. CALCULATION OF UNOCCUPIED ORBITALS
The method reduces basis-function requirements for a target precision and maintains high parallel efficiency, while large systems become dominated by cubic linear-algebra operations.
- Performance results: The wavelet method reduces the number of degrees of freedom needed to attain a given absolute precision compared with a plane-wave code.This reduction lowers memory requirements and floating-point operations for molecular calculations.
- Performance results: 88% overall parallel efficiency is maintained even for large systems using many processors.
- Performance results: For relatively small systems, the Poisson solver dominates execution time, whereas linear-algebra operations dominate for large systems.The latter include matrix products and decompositions associated with orthogonality constraints and orthogonalization.
- Performance results: The dominant large-system linear-algebra operations scale cubically with the number of atoms.The overlap-matrix calculation and orbital orthogonality transformation produce this scaling, with degrees of freedom typically exceeding the number of orbitals.
- Method and implementation: The method analytically computes matrix elements, kinetic energy, and nonlocal pseudopotential operators, while other operations use short-range-filter convolutions.The implementation is integrated into ABINIT and distributed under the GNU-GPL license.