PaperlayerView exact revision on arXiv

Source-linked AI summary

Higher-order adaptive finite-element methods for Kohn-Sham density functional theory

Phani Motamarri, Michael R Nowak, Kenneth Leiter, Jaroslaw Knap, Vikram Gavini

arXiv:1207.0167v5physics.comp-ph

TL;DR

Kohn-Sham DFT needs efficient, accurate real-space discretizations that avoid the basis-size limitations of linear finite-elements. The paper combines adaptive higher-order spectral finite-elements with efficient eigenvalue solvers, achieving roughly 1000-fold savings over linear elements, competitiveness with plane waves for non-periodic pseudopotential calculations, and calculations on a 1,688-atom cluster.

  • Problem

    Linear finite-elements require many basis functions to reach chemical accuracy, motivating more efficient higher-order real-space discretizations for Kohn-Sham DFT.

  • Method

    The paper combines a-priori mesh adaptation with spectral finite-elements, Gauss-Lobatto quadrature, and Chebyshev acceleration for the occupied eigenspace.

  • Results

    1000-fold computational savings over linear finite-elements are observed at chemical accuracy, while sixth-order elements compete with plane waves for non-periodic pseudopotential calculations and reach 1,688 atoms.

  • Takeaways & Limitations

    Higher-order finite-elements can provide chemically accurate Kohn-Sham calculations with substantially lower cost than linear elements and practical capability for large non-periodic systems.

  • Takeaways & Limitations

    Using the same finite-element discretization for wavefunctions and electrostatic potential makes discretization convergence non-variational in general, and all-electron finite-elements are about ten times slower than optimized Gaussian bases in the reported benchmarks.

Abstract

from arXiv · show

We present an efficient computational approach to perform real-space electronic structure calculations using an adaptive higher-order finite-element discretization of Kohn-Sham density-functional theory (DFT). To this end, we develop an a-priori mesh adaption technique to construct a close to optimal finite-element discretization of the problem. We further propose an efficient solution strategy for solving the discrete eigenvalue problem by using spectral finite-elements in conjunction with Gauss-Lobatto quadrature, and a Chebyshev acceleration technique for computing the occupied eigenspace. The proposed approach has been observed to provide a staggering 100-200 fold computational advantage over the solution of a generalized eigenvalue problem. Using the proposed solution procedure, we investigate the computational efficiency afforded by higher-order finite-element discretization of the Kohn-Sham DFT problem. Our studies suggest that staggering computational savings of the order of 1000 fold relative to linear finite-elements can be realized, for both all-electron and local pseudopotential calculations. On all the benchmark systems studied, we observe diminishing returns in computational savings beyond the sixth-order for accuracies commensurate with chemical accuracy. A comparative study of the computational efficiency of the proposed higher-order finite-element discretizations suggests that the performance of finite-element basis is competing with the plane-wave discretization for non-periodic local pseudopotential calculations, and compares to the Gaussian basis for all-electron calculations within an order of magnitude. Further, we demonstrate the capability of the proposed approach to compute the electronic structure of a metallic system containing 1688 atoms using modest computational resources, and good scalability of the present implementation up to 192 processors.

1 Introduction

The paper develops adaptive higher-order finite-element methods for efficient real-space Kohn-Sham DFT, addressing the large basis requirements of linear elements while retaining finite-elements’ geometric flexibility. Its benchmarks report substantial savings, competitive performance against plane waves for non-periodic pseudopotential systems, and practical scaling to large clusters.

  • Motivation: Higher-order finite-elements address the chemical-accuracy inefficiency of linear elements, which can require about 100,000 basis functions per atom.Finite-elements remain attractive because they support unstructured coarse-graining, complex geometries, boundary conditions, and parallel scalability.
  • Methods: The approach combines a-priori mesh adaptation with spectral finite-elements, Gauss-Lobatto quadrature, and Chebyshev acceleration for the occupied eigenspace.The mesh distribution is derived by minimizing estimated ground-state discretization error subject to a fixed element count, while the solver strategy addresses the generalized eigenproblem.
  • Benchmarks: The study evaluates convergence rates on all-electron boron and methane systems and local-pseudopotential barium and periodic calcium benchmarks.The analysis addresses higher-order finite-element discretization for Kohn-Sham DFT involving Coulomb-singular potentials.
  • Computational efficiency: 1000-fold computational savings over linear finite-elements are observed for both all-electron and local pseudopotential calculations at chemical accuracy, with diminishing returns beyond sixth order.The diminishing returns arise because increasing per-basis-function costs offset the degree-of-freedom advantage beyond sixth order.
  • Comparative benchmarks: Higher-order finite-elements use less CPU time than plane waves for comparable ground-state accuracy in non-periodic aluminium clusters, and compute a 1,688-atom cluster that ABINIT could not simulate because of memory requirements.The comparison uses sixth-order spectral elements and local pseudopotentials.
  • Comparative benchmarks: For all-electron benchmarks, finite-elements require about ten times more solution time than Gaussian basis calculations, partly because optimized Gaussian functions need fewer basis functions.The authors note that improved parallel scalability and enriched finite-element shape functions could reduce this gap.

2 Formulation

The formulation expresses Kohn-Sham DFT as a variational ground-state problem and a self-consistent nonlinear eigenvalue problem, then discretizes the fields with finite elements. Electrostatic interactions are localized through an electrostatic potential, while the discrete problem is generally a generalized Hermitian eigenproblem whose convergence can be non-variational when wavefunctions and potentials share the same discretization.

  • 2.1 Kohn-Sham variational problem: The Kohn-Sham energy functional combines noninteracting kinetic, exchange-correlation, Hartree, external-potential, and nuclear-repulsion energies.The study uses the local-density approximation for exchange-correlation energy and notes extensions to LSDA and GGA.
  • 2.1 Kohn-Sham variational problem: The electrostatic energy terms depend on whether the system is finite or periodic, with integrations and nuclear sums taken over the corresponding real-space domains.Finite systems use integrals over R3, whereas periodic calculations integrate over the unit cell with the appropriate extended electrostatic contribution.
  • 2.1 Kohn-Sham variational problem: Electrostatic interactions are reformulated as a local variational problem in the total electrostatic potential, replacing nonlocal real-space terms with a Poisson-based formulation.The potential includes contributions from electron density and nuclear charge distributions, with nuclear self-energy treated as an inconsequential constant that is subtracted numerically.
  • 2.1 Kohn-Sham variational problem: The ground-state energy and electron density are obtained by minimizing the energy over orthonormal single-electron wavefunctions and suitable function spaces.Non-periodic and periodic problems use different bounded-domain spaces, including H1_0(Ω) and periodic H1 spaces.
  • 2.2 Kohn-Sham eigenvalue problem: Stationarity of the variational problem yields a nonlinear self-consistent Kohn-Sham eigenvalue problem whose effective potential depends on electron density.Periodic systems additionally exploit lattice periodicity and Bloch solutions, while all systems require appropriate boundary conditions.
  • 2.3 Discrete Kohn-Sham eigenvalue problem: Finite-element expansions approximate both electronic wavefunctions and electrostatic potential in a finite-dimensional subspace, producing discrete Hamiltonian and overlap matrices.The resulting problem is a generalized eigenvalue problem because the finite-element basis is non-orthogonal.
  • 2.3 Discrete Kohn-Sham eigenvalue problem: Using the same finite-element discretization for wavefunctions and electrostatic potential makes discretization convergence non-variational in general.The drawback can be mitigated by refining the electrostatic-field mesh or using higher-order polynomials for the potential than for the wavefunctions.

3 A-priori mesh adaption

The paper develops an a-priori mesh-adaptation scheme for Kohn-Sham DFT that minimizes finite-element energy error under a fixed element budget. The resulting mesh distribution is derived from energy-error bounds expressed through wavefunction and electrostatic-potential approximation errors.

  • 3 A-priori mesh adaption: The proposed a-priori scheme minimizes Kohn-Sham finite-element approximation error for a fixed number of mesh elements.The method derives an energy-error bound and uses it to determine the mesh-size distribution.
  • 3.1 Estimate of energy error: The ground-state energy error is bounded in terms of approximation errors in the electronic wavefunctions and electrostatic potential.The analysis expands the discrete energy around the continuous solution and relates the resulting error to these electronic fields.
  • 3.2 Optimal mesh-size distribution: The optimal mesh distribution is obtained by minimizing the error estimate subject to a fixed total number of elements.A variational problem, its Euler-Lagrange equation, and the resulting distribution determine the target element size h(r).
  • 3.2 Optimal mesh-size distribution: The mesh construction uses asymptotic knowledge of the wavefunctions and electrostatic potential to apply the derived coarse-graining rate.This provides the spatial information needed to construct meshes for the subsequent problem classes.

4 Numerical implementation

The implementation combines higher-order and spectral finite elements with Gauss-Lobatto quadrature to simplify the Kohn-Sham eigenproblem, then uses eigenspace methods within self-consistent iterations. The resulting strategy compares standard eigensolvers with Chebyshev filtering and addresses charge-sloshing through finite-temperature occupations.

  • 4.1 Higher-order finite-element discretizations: The implementation uses C0 finite elements from linear tetrahedra through degree-eight spectral hexahedra to target chemical accuracy more efficiently.The studied elements include TET4, HEX27, HEX125SPECT, HEX343SPECT, and HEX729SPECT.
  • 4.1 Higher-order finite-element discretizations: Gauss-Lobatto-Legendre quadrature makes the spectral-element overlap matrix diagonal, enabling a trivial transformation from the generalized to the standard eigenproblem.Quadrature points coincide with nodal points, causing off-diagonal mass-matrix entries to vanish.
  • 4.2 Numerical solution strategy: The occupied eigenspace is computed either by explicit eigenvector solvers or by Chebyshev filtering during self-consistent-field iterations.The explicit route compares generalized and transformed standard Hermitian eigenproblems, while the filtering route accelerates subspace computation.
  • 4.2.1.1 Explicit computation of eigenvectors: Transforming the Kohn-Sham problem to a standard eigenproblem yields a 10-fold speedup for graphene and a 25-fold speedup for the aluminum cluster.The Krylov-Schur method uses the lowest-eigenvalue vector from the preceding self-consistent iteration as its input approximation.
  • 4.2.2 Finite temperature smearing: Finite-temperature Fermi-Dirac occupations are used to mitigate charge sloshing in systems with very small band gaps or degenerate Fermi-level states.Charge sloshing can produce large changes in spatial charge density between self-consistent iterations.

5.1 Rates of convergence

The convergence study evaluates higher-order finite elements across all-electron, local-pseudopotential, cluster, and periodic-crystal benchmarks. Higher-order elements achieve near-optimal convergence and substantially reduce the element counts needed for chemical accuracy compared with linear elements.

  • Benchmark systems: The study compares convergence for boron and methane all-electron systems with local-pseudopotential barium and periodic calcium benchmarks.The benchmark set spans isolated atoms and molecules, a non-periodic cluster, and a bulk FCC crystal.
  • Convergence behavior: Close to optimal convergence rates O(h2k) are obtained across the studied finite-element discretizations, despite nonlinear equations and Coulomb-singular nuclear potentials.Here, k is the polynomial degree.
  • Boron atom: Five million linear TET4 elements are needed for boron at relative error 10^-2, whereas a few hundred HEX125SPECT or HEX343SPECT elements reach 10^-4.Higher accuracies require only a few thousand higher-order elements.
  • Barium cluster: For the barium cluster, tetrahedral and hexahedral elements show near-optimal convergence, while its finite-element reference energy is −0.6386307998 Ha versus −0.638627743 Ha from ABINIT.The comparison uses the stated plane-wave calculation with a 30 Ha cutoff and 80 a.u. cell size.

5.2 Computational cost

Higher-order finite elements substantially reduce computational cost across benchmark Kohn-Sham DFT calculations, with benefits increasing toward chemical accuracy and diminishing beyond sixth order. The approach also remains competitive with plane waves for non-periodic local-pseudopotential calculations and supports a 1688-atom aluminum system.

  • 5.2.1 Benchmark systems: A thousand-fold computational advantage over linear TET4 elements is obtained even at relative errors of 10^-2.
  • 5.2.1 Benchmark systems: At relative error 10^-5, HEX125SPECT is nearly ten times faster than HEX27, while HEX343SPECT is more than three times faster than HEX125SPECT.
  • 5.2.1 Benchmark systems: At chemical accuracy, computational efficiency improves significantly through sixth-order elements, with diminishing returns beyond sixth order.Increasing Hamiltonian-matrix construction costs increasingly dominate total time at higher orders.
  • 5.2.2 Larger systems: Higher-order finite elements achieve comparable, and in one case two-fold better, computational efficiency than plane waves for non-periodic local-pseudopotential calculations.
  • 5.2.2 Larger systems: A finite-element calculation was performed for a 1688-atom aluminum cluster, while the corresponding plane-wave calculation was beyond available computational resources.
  • 5.2.2 Larger systems: Gaussian calculations were faster in the all-electron benchmarks, partly because their basis sets were highly optimized for the specific material systems.

5.3 Scalability of finite-element basis:

The implementation exhibits nearly linear strong scaling on a fixed aluminum-cluster problem as processor count increases. With 96 times more processors, it reaches 91.4% efficiency.

  • 5.3 Scalability of finite-element basis:: The implementation achieves an 87.82-fold speedup and 91.4% efficiency for a 96-fold increase in processors.The fixed problem is an aluminum 3 × 3 × 3 cluster with 3.91 million degrees of freedom; timings are referenced to two processors.
  • 5.3 Scalability of finite-element basis:: Strong scaling is measured on a fixed-size aluminum cluster using HEX125SPECT elements with 3.91 million degrees of freedom.
  • 5.3 Scalability of finite-element basis:: The observed speedup is described as almost linear over the tested processor range.

6 Conclusions

The study demonstrates that adaptive higher-order spectral finite elements can substantially improve Kohn-Sham DFT efficiency while retaining accurate discretization and broad applicability.

  • 6 Conclusions: The transformed standard eigenvalue problem and Chebyshev filtering provide the computational framework for evaluating higher-order finite-element efficiency.The approach combines adaptive meshes, spectral elements, Gauss-Lobatto quadrature, and Chebyshev filtering.
  • 6 Conclusions: Optimal convergence rates are obtained across all tested finite-element orders, including all-electron calculations with Coulomb-singular potentials.The meshes are constructed using error estimates and far-field asymptotics for benchmark systems.
  • 6 Conclusions: 1000-fold CPU-time savings are achieved with sixth-order hexahedral spectral elements over linear tetrahedral elements.The comparison covers benchmark all-electron and local-pseudopotential calculations.
  • 6 Conclusions: Finite elements compete with plane waves for non-periodic local-pseudopotential calculations, while all-electron calculations are within an order of magnitude of Gaussian bases.The study also computes a 1,688-atom aluminium cluster that was not feasible with ABINIT because of memory requirements.
  • 6 Conclusions: The study identifies higher-order spectral finite elements as a promising direction, with further reductions possible through single-atom enrichment and future linear-scaling formulations.These extensions are presented as future work rather than demonstrated results.

A Discrete formulation of electrostatic interactions in all-electron calculations

The discrete electrostatic formulation represents electron, nuclear, and self-energy contributions in the finite-element basis and yields convergent total electrostatic energies despite divergent individual self-energy terms.

  • A Discrete formulation of electrostatic interactions in all-electron calculations: The finite-element formulation computes electrostatic potentials from the electron density and nuclear charge distribution, with nuclear point charges located on mesh nodes.The mesh discreteness regularizes the nuclear potential representation.
  • A Discrete formulation of electrostatic interactions in all-electron calculations: The electrostatic energy combines electron-electron, electron-nuclear, nuclear-nuclear, and nuclear self-energy contributions.The formulation separates the additional nuclear interaction and self-energy terms into contributions labeled (a) and (b).
  • A Discrete formulation of electrostatic interactions in all-electron calculations: Despite divergent terms (a) and (b) under mesh refinement, the total electrostatic energy converges because the divergent self-energy contributions cancel.Methane convergence is reported as close to optimal, with extrapolated E0 = −23.79671760794.

B Trade-offs in higher-order elements: Source of diminishing returns with increasing order

For chemical-accuracy-level errors, increasing element order improves efficiency through sixth order, after which rising matrix-construction costs erase the savings from fewer basis functions.

  • B Trade-offs in higher-order elements: Source of diminishing returns with increasing order: Beyond sixth order, higher-order elements provide diminishing computational returns for the benchmark barium-cluster problem.The comparison uses meshes with similar relative ground-state-energy errors of about 10^-5.
  • B Trade-offs in higher-order elements: Source of diminishing returns with increasing order: Increasing polynomial degree reduces basis-function counts but does not make octic elements cheaper than hexic elements because Hamiltonian construction becomes more expensive.The cost rises with the number of basis functions per element and the quadrature order.
  • B Trade-offs in higher-order elements: Source of diminishing returns with increasing order: For the studied systems, Hamiltonian construction dominates increasing-order costs; in larger systems, Chebyshev-vector orthogonalization could shift the diminishing-return point beyond sixth order.The latter possibility is explicitly conditional on orthogonalization becoming the dominant SCF cost.

C Accuracy of Gauss-Lobatto-Legendre quadrature

The study evaluates Gauss-Lobatto-Legendre quadrature for overlap matrices in spectral finite elements and examines its convergence and comparison with standard Gauss quadrature.

  • C Accuracy of Gauss-Lobatto-Legendre quadrature: GLL quadrature produces a diagonal overlap matrix, converting the finite-element generalized eigenvalue problem into a standard eigenvalue problem suitable for Chebyshev filtering.This is the computational rationale for using reduced-order quadrature in the spectral-element formulation.
  • C Accuracy of Gauss-Lobatto-Legendre quadrature: Hydrogen-atom ground-state energies converge under GLL quadrature across a sequence of refined spectral-element meshes.The test uses HEX125SPECTRAL elements on a spherical domain of radius 20 a.u.
  • C Accuracy of Gauss-Lobatto-Legendre quadrature: GLL and Gauss quadrature are compared for methane and a 2×2×2 barium cluster using coarse and relatively fine meshes across different discretization orders.The results are tabulated for both all-electron and local-pseudopotential benchmark systems.
Loading 1207.0167v5…