Source-linked AI summary
Adaptive local basis set for Kohn-Sham density functional theory in a discontinuous Galerkin framework I: Total energy calculation
Lin Lin, Jianfeng Lu, Lexing Ying, E. Weinan
TL;DR
Uniform Kohn-Sham discretizations require many basis functions near nuclei, while atomic-orbital approaches need fine tuning for high accuracy. The paper introduces adaptive local DG basis functions that incorporate orbital oscillations and environmental effects, achieving less than 1meV accuracy with 4–40 basis functions per atom and supporting systems with thousands of atoms.
Problem
Uniform discretization requires many basis functions per atom, while atomic-orbital approaches generally require fine tuning to reach high accuracy.
Method
The method adaptively constructs localized, discontinuous basis functions from local environments and uses the DG framework to reconstruct continuous Kohn-Sham orbitals and electron density.
Results
Less than 1meV total-energy accuracy is achieved with 4–40 basis functions per atom, while the parallel implementation handles systems with thousands of atoms.
Takeaways & Limitations
Environmental effects can be built systematically into localized basis functions while retaining high-accuracy total-energy and electron-density calculations.
Takeaways & Limitations
Force calculations and geometry-based simulations require access to adaptive basis functions and their atomic-position derivatives, including Pulay forces.
Abstract
from arXiv · showhide
Kohn-Sham density functional theory is one of the most widely used electronic structure theories. In the pseudopotential framework, uniform discretization of the Kohn-Sham Hamiltonian generally results in a large number of basis functions per atom in order to resolve the rapid oscillations of the Kohn-Sham orbitals around the nuclei. Previous attempts to reduce the number of basis functions per atom include the usage of atomic orbitals and similar objects, but the atomic orbitals generally require fine tuning in order to reach high accuracy. We present a novel discretization scheme that adaptively and systematically builds the rapid oscillations of the Kohn-Sham orbitals around the nuclei as well as environmental effects into the basis functions. The resulting basis functions are localized in the real space, and are discontinuous in the global domain. The continuous Kohn-Sham orbitals and the electron density are evaluated from the discontinuous basis functions using the discontinuous Galerkin (DG) framework. Our method is implemented in parallel and the current implementation is able to handle systems with at least thousands of atoms. Numerical examples indicate that our method can reach very high accuracy (less than 1meV) with a very small number ($4\sim 40$) of basis functions per atom.
1. Introduction
The paper targets the large basis requirements of uniform Kohn-Sham discretizations and the tuning limitations of atomic-orbital approaches. It proposes adaptive, environmentally informed discontinuous Galerkin basis functions that achieve high accuracy with few basis functions per atom.
- Motivation: 500–5000 basis functions per atom are typically required by standard planewave discretization at 15Ry–90Ry kinetic-energy cutoffs.Other uniform discretization methods generally require even more basis functions per atom.
- Motivation: Rapid orbital oscillations near nuclei drive the large basis requirement, while atomic orbitals capture a significant part of these oscillations.The orbitals are smoother in the interstitial region, motivating localized basis construction.
- Limitations of prior approaches: Atomic-orbital and mixed-basis methods still require parameter fine tuning because their basis functions omit environmental effects.In some circumstances, these methods also retain a large number of basis functions per atom.
- Proposed approach: The proposed method adaptively constructs discontinuous local basis functions from atomic configurations and uses the DG framework to approximate continuous orbitals.The basis functions incorporate environmental effects while remaining localized in real space.
- Reported capability: Less than 1meV total-energy accuracy is obtained with 4–40 basis functions per atom, and the parallel implementation handles systems with thousands of atoms.The paper organizes the method around DG formulation, adaptive basis construction, implementation, and performance evaluation.
2. Discontinuous Galerkin framework for Kohn-Sham density functional theory
The DG formulation discretizes the nonlinear Kohn-Sham problem using element-local, discontinuous basis functions and an interior-penalty treatment of interfaces. Self-consistent iteration updates the electron density through generalized eigenvalue solves.
- KSDFT formulation: Kohn-Sham DFT seeks the total energy and self-consistent electron density for a specified atomic configuration.The effective potential depends on the density, making the Kohn-Sham equation a nonlinear eigenvalue problem.
- Self-consistent iteration: The self-consistent iteration repeatedly constructs an effective Hamiltonian, evaluates the density from its eigenfunctions, and mixes densities until convergence.The implementation uses Anderson mixing, although linear mixing is presented as the simplest choice.
- Computational focus: The electron-density evaluation is the main bottleneck in the self-consistent iteration and is the focus of the paper’s numerical discretization.The method targets efficient and accurate density evaluation for a given effective potential.
- DG framework: DG methods permit discontinuous basis functions across element interfaces, providing flexibility in constructing the approximation space.The paper uses the interior-penalty method, whose penalty parameter controls jumps across element surfaces and guarantees stability.
- Domain and partition: The computational domain uses periodic boundary conditions and is partitioned into quasi-uniform rectangular elements with corresponding surfaces.Elements are typically chosen with equal size, such as integer multiples of a crystal’s conventional cell.
- Discrete eigenproblem: Element-local basis functions are extended by zero outside their elements, assembled into a global approximation space, and used to derive a generalized eigenvalue problem.The resulting matrices are the DG stiffness matrix and mass matrix; an identity mass matrix yields a standard eigenvalue problem.
3. Basis functions adapted to the local environment
The method constructs localized, discontinuous basis functions from eigenfunctions of locally extended Hamiltonians, incorporating both nuclear oscillations and environmental effects. SVD filtering orthonormalizes the basis, while completeness as the local basis size grows enables systematic accuracy improvement.
- Motivation: Uniform fine discretization is costly because Kohn-Sham orbitals oscillate near nuclei but are smooth in interstitial regions.This motivates adapting the basis to local structure rather than resolving the entire domain uniformly.
- Motivation: Unlike isolated-atom approaches, the proposed basis functions incorporate information from each element’s local environment and remain systematically improvable.The DG framework permits these locally adapted functions to form a complete global basis.
- Adaptive local basis construction: For each element, eigenfunctions of a periodically restricted Hamiltonian on an extended region are restricted to the element and zero-extended globally.The extended region includes a buffer area intended to capture nearby environmental effects and local orbital behavior.
- Adaptive local basis construction: Periodic boundary conditions on each extended region capture bulk behavior and facilitate use of existing DFT algorithms, although optimal boundary conditions remain open.Other choices, including Neumann conditions, are possible.
- Accuracy and conditioning: As the number of retained local functions increases, the restricted functions become complete on each element, allowing systematic accuracy improvement.This follows from completeness of the extended-region Hamiltonian eigenfunctions.
- Accuracy and conditioning: SVD filtering removes numerically degenerate directions and produces an orthonormal basis, converting the generalized eigenvalue problem into a standard one.The filtering also controls conditioning as the candidate basis size increases.
- Algorithm: The overall algorithm alternates adaptive local-basis construction, DG solution of the Kohn-Sham problem, and electron-density reconstruction.The basis is rebuilt from the effective Hamiltonian associated with the current trial density.
4. Implementation details
The implementation uses distinct global, extended-element, and element grids, with interpolation and local numerical integration connecting them. Parallel communication is localized, but the current eigensolver treats the sparse DG stiffness matrix as dense and requires redistribution.
- 4.1. Grids and interpolation: The implementation distinguishes global, extended-element, and element domains, each using grids suited to its quantities and operations.The global and extended-element grids are uniform Cartesian grids, while elements use three-dimensional LGL grids.
- 4.1. Grids and interpolation: Interpolation transfers density and potentials from the global grid, basis functions from extended elements to elements, and assembled orbitals back to the global grid.Fourier interpolation is used for extended-element-to-element transfer, while Lagrange interpolation and averaging are used for element-to-global transfer.
- 4.1. Grids and interpolation: Non-local pseudopotentials are generated directly in real space on both extended-element and element grids, avoiding further interpolation between those grids.They enter both local Hamiltonian solves and elementwise numerical integration.
- 4.2. DG implementation: Planewaves discretize each extended-element Hamiltonian, LOBPCG diagonalizes it, and the resulting eigenfunctions are restricted to elements before SVD filtering.The filtering uses a discrete weighted ℓ2 inner product with LGL weights and yields a standard eigenvalue problem.
- 4.2. DG implementation: The DG stiffness matrix combines element-local volume terms, face integrals, and localized pseudopotential terms, producing block-diagonal and sparse components.Volume contributions are nonzero only within individual elements, while pseudopotential contributions require overlapping supports.
- 4.3. Parallel implementation: The current implementation does not exploit stiffness-matrix sparsity: ScaLAPACK’s dense eigensolver requires storing the matrix densely and redistributing it.The authors plan to replace this with solvers that leverage sparsity.
- 4.3. Parallel implementation: The algorithm is fully parallel, with processors computing local Hamiltonians independently after exchanging pseudopotentials from nearby elements.DG assembly requires neighboring basis functions and basis functions overlapping nearby non-local pseudopotentials.
5. Numerical examples
Numerical tests on Na and Si assess accuracy, dimensionality, penalty-parameter effects, and computational efficiency of adaptive local basis functions in the DG framework.
- Overall performance: 10^-6 au per atom accuracy is achieved in ground-state calculations using a small number of adaptive local basis functions.The study also reports significant speedup for a 128-atom Na system and implementation on systems with thousands of atoms.
- Test setup: The tests use Na and Si as representative metallic and insulating systems, respectively, with LDA and HGH pseudopotentials.Calculations use atomic units and are performed on the Hopper system at NERSC.
- Disordered quasi-1D systems: 4.3 × 10^-7 au or 0.01 meV total-energy error per atom is obtained for quasi-1D Na with 10 basis functions per atom and buffer size 1.00.With buffer size 0.50, the error is already below 10^-3 au using 4 basis functions per atom.
- Disordered quasi-1D systems: 7.8×10^-8 au or 0.002 meV total-energy error per atom is obtained for quasi-1D Si with 8 basis functions per atom and buffer size 1.00.The reported 20 ∼40 basis functions per atom for Gaussian or numerical atomic orbitals contrasts with the tested adaptive basis sizes.
- Disordered quasi-2D and 3D systems: 1.2 × 10^-3 au error with 24 basis functions per atom at buffer size 0.50 decreases to 5.6 × 10^-6 au or 0.15 meV with 42 functions and buffer size 1.00 in bulk 3D Na.For quasi-2D Na, the error reaches 2.8 × 10^-6 au with 16 basis functions per atom and buffer size 1.00.
- Penalty parameter: The penalty parameter affects both stability and error: errors scale empirically as α0.66 for Na and α0.58 for Si up to α = 640.The DG formulation can become unstable below a problem-dependent threshold, including α = 5 for Na and α = 10 for Si.
- Computational efficiency: Adaptive basis dimension reduction yields significant computational-time savings, while the reported comparison excludes the cost of solving for the adaptive local basis functions.The efficiency comparison contrasts global-domain LOBPCG costs with parallel DG eigensolver wall-clock times.
6. Discussion and Conclusion
The paper demonstrates adaptive local basis functions within a discontinuous Galerkin framework for accurate total-energy calculations with few basis functions per atom. The implementation handles systems of thousands of atoms, while force-related extensions and further performance improvements remain ongoing.
- Adaptive local basis functions efficiently calculate total energy and electron density with high accuracy using very few basis functions per atom.The discontinuous basis functions reconstruct continuous Kohn-Sham orbitals and electron density through the DG framework, while incorporating environmental effects automatically.
- The current implementation performs total-energy calculations for systems consisting of thousands of atoms.
- Force calculations, geometry optimization, and ab initio molecular dynamics require access to adaptive basis functions and their atomic-position derivatives.The paper identifies control of the resulting Pulay force through optimized local basis functions as ongoing work.