PaperlayerView exact revision on arXiv

Source-linked AI summary

DFT-FE -- A massively parallel adaptive finite-element code for large-scale density functional theory calculations

Phani Motamarri, Sambit Das, Shiva Rudraraju, Krishnendu Ghosh, Denis Davydov, Vikram Gavini

arXiv:1903.10959v3physics.comp-ph

TL;DR

Large-scale Kohn–Sham DFT calculations require substantial computational resources, motivating scalable methods for systems with many thousands of atoms. DFT-FE combines a local real-space adaptive spectral finite-element formulation with scalable eigensolver strategies, achieving accurate results and substantially better efficiency than QE for larger systems.

  • Problem

    Large-scale DFT calculations are computationally demanding, while Kohn–Sham eigensolver costs typically scale cubically with system size.

  • Method

    DFT-FE uses a local real-space adaptive spectral finite-element discretization with Chebyshev filtering, Cholesky-based orthonormalization, spectrum splitting, and mixed precision.

  • Results

    DFT-FE agrees excellently with established DFT codes and is 4.5–12× more efficient than QE in CPU-time for systems containing 10,659–20,470 electrons.

  • Takeaways & Limitations

    DFT-FE provides accurate, massively parallel DFT calculations for pseudopotential and all-electron systems with periodic, non-periodic, or semi-periodic boundaries.

  • Takeaways & Limitations

    Octree-based non-conforming refinement can only approximately realize target mesh sizes and may create excess degrees of freedom, especially for all-electron calculations.

Abstract

from arXiv · show

We present an accurate, efficient and massively parallel finite-element code, DFT-FE, for large-scale ab-initio calculations (reaching $\sim 100,000$ electrons) using Kohn-Sham density functional theory (DFT). DFT-FE is based on a local real-space variational formulation of the Kohn-Sham DFT energy functional that is discretized using a higher-order adaptive spectral finite-element (FE) basis, and treats pseudopotential and all-electron calculations in the same framework, while accommodating non-periodic, semi-periodic and periodic boundary conditions. We discuss the main aspects of the code, which include, the various strategies of adaptive FE basis generation, and the different approaches employed in the numerical implementation of the solution of the discrete Kohn-Sham problem that are focused on significantly reducing the floating point operations, communication costs and latency. We demonstrate the accuracy of DFT-FE by comparing the energies, ionic forces and periodic cell stresses on a wide range of problems with popularly used DFT codes. Further, we demonstrate that DFT-FE significantly outperforms widely used plane-wave codes---both in CPU-times and wall-times, and on both non-periodic and periodic systems---at systems sizes beyond a few thousand electrons, with over $5-10$ fold speedups in systems with more than 10,000 electrons. The benchmark studies also highlight the excellent parallel scalability of DFT-FE, with strong scaling demonstrated on up to 192,000 MPI tasks.

1. Introduction

DFT-FE addresses the computational and scalability limits of large-scale Kohn-Sham DFT through an adaptive real-space finite-element formulation and optimized eigensolver implementation. The paper validates its accuracy and demonstrates efficient, massively parallel calculations across diverse systems.

  • Motivation: Large-scale DFT is needed for materials modeling, but computational costs routinely limit accurate calculations to a few thousand electrons.This restricts simulations involving many atoms, long molecular-dynamics trajectories, and numerous relaxation steps.
  • Motivation: Plane-wave bases restrict domains to periodic systems and lack good parallel scalability, while atomic-orbital bases are not systematically convergent for generic materials.
  • Motivation: Cubic-scaling diagonalization of the Kohn-Sham Hamiltonian creates another major computational barrier as electron number increases.
  • Approach: DFT-FE combines adaptive spectral finite elements with Chebyshev filtering, Cholesky-based orthonormalization, spectrum splitting, and mixed precision to reduce computational costs.The implementation targets floating-point operations, communication costs, and latency.
  • Approach: The code uses a local real-space variational formulation that unifies pseudopotential and all-electron calculations and supports periodic, non-periodic, and semi-periodic boundary conditions.The formulation also supports configurational-force calculations for ionic forces and stresses.
  • Performance: Domain decomposition combined with moderate band parallelization achieves about 50% scaling efficiency on 102,400 MPI tasks for a system with about 40,000 electrons.The reported scalability is attributed to the locality of the finite-element basis and communication-reducing implementations.
  • Validation: Validation compares DFT-FE with established plane-wave and all-electron codes using energies, ionic forces, stresses, relaxations, molecular dynamics, and band structures.

2. Real space Kohn-Sham DFT formulation

The real-space formulation expresses Kohn-Sham DFT through local variational and Poisson problems while supporting pseudopotential and all-electron treatments under multiple boundary conditions. It also provides a unified variational route to ionic forces and stresses.

  • 2.1. Governing equations in DFT: Kohn-Sham DFT determines ground-state properties by solving for the lowest eigenvalues and eigenfunctions of a nonlinear effective Hamiltonian.
  • 2.1. Governing equations in DFT: Electron density is constructed from canonical wavefunctions and orbital occupancies, with the Fermi energy fixed by the total electron number.DFT-FE uses a Fermi–Dirac occupancy function to regularize systems with many states near the Fermi energy.
  • 2.1. Governing equations in DFT: The effective single-electron potential combines exchange-correlation and electrostatic contributions.The exchange-correlation potential is the variational derivative of the exchange-correlation energy; LDA and GGA forms are incorporated.
  • 2.1. Governing equations in DFT: Electrostatic energies include electron–electron, electron–nucleus, and nucleus–nucleus interactions, with their definitions adapted to periodic or non-periodic domains.
  • 2.1. Governing equations in DFT: The pseudopotential operator separates local and nonlocal terms and uses norm-conserving projector representations, including Troullier–Martins and ONCV forms.
  • 2.1. Governing equations in DFT: Extended electrostatic interactions are reformulated as local variational problems by computing associated potentials through Poisson equations.
  • 2.1. Governing equations in DFT: All-electron calculations arise by removing the pseudopotential terms, while the same equations apply to periodic and non-periodic systems with appropriate boundary conditions.
  • 2.3. Configurational Forces: Configurational forces are derivatives of the Kohn-Sham energy with respect to material-point positions, providing ionic forces and stresses while inherently including Pulay corrections.

3. Numerical implementation

DFT-FE is implemented on the deal.II finite-element library, using its finite-element constructs, adaptive-mesh architecture, and parallel vector objects.

  • Implementation: DFT-FE uses deal.II for finite-element constructs, adaptive mesh refinement, and parallel vector operations.

3.1. Adaptive mesh refinement

DFT-FE supports octree-based adaptive finite-element meshes, using user-defined or automatic strategies to concentrate resolution where needed. The automatic strategy avoids repeated Kohn–Sham solves during refinement and achieves substantially lower degrees of freedom while retaining systematic energy and force convergence.

  • Adaptive mesh construction: Octree-based hexahedral meshes provide adaptive spatial resolution for DFT-FE calculations.The mesh generator uses the p4est library through deal.II.
  • User-defined adaptive mesh refinement: User-defined refinement marks cells near atoms within radius ratom and refines them from hbase toward hatom.The resulting mesh coarsens away from atomic neighborhoods; hbase, hatom, and ratom are user inputs.
  • User-defined adaptive mesh refinement: The target mesh distribution h(x) is estimated by minimizing energy approximation error subject to a fixed number of elements.Numerically constructed single-atom wavefunctions provide practical estimates for hbase, hatom, and hfine.
  • Automatic adaptive mesh refinement: AAMR constructs the mesh before SCF using single-atom Kohn–Sham wavefunctions, avoiding expensive repeated large-scale Kohn–Sham solves during refinement.Its error indicator is based on an energy estimate involving discrete wavefunctions.
  • Automatic adaptive mesh refinement: 20 times fewer degrees of freedom were obtained with AAMR than UDAMR for all-electron SiF4 at ground-state energy and force errors of O(10^-4).For periodic all-electron Si diamond, AAMR used 3 times fewer degrees of freedom at ground-state energy errors of O(10^-4).
  • Automatic adaptive mesh refinement: AAMR attained force discretization errors as low as O(10^-5) Ha/Bohr, although errors stagnated on increasingly finer meshes.The stagnation is attributed to the a priori error indicator using single-atom wavefunctions.

3.2. SCF Algorithm

DFT-FE solves the discrete nonlinear Kohn–Sham problem self-consistently using adaptive higher-order spectral finite elements and Chebyshev-filtered subspace iteration. Its computational design targets the different scaling regimes of filtering, orthonormalization, and Rayleigh–Ritz operations.

  • SCF solution: The SCF procedure solves the discrete nonlinear Hermitian Kohn–Sham eigenproblem together with Poisson equations.Adaptive higher-order spectral finite elements are combined with ChFSI to evaluate the occupied eigenspace.
  • Computational scaling: Chebyshev filtering scales quadratically with atom number, whereas Cholesky-based orthonormalization and Rayleigh–Ritz scale cubically.Filtering dominates small-to-medium systems, while orthonormalization and Rayleigh–Ritz dominate larger systems.

3.3. Chebyshev filtering

DFT-FE uses Chebyshev filtering to approximate the occupied eigenspace within each self-consistent iteration, followed by orthonormalization and Rayleigh–Ritz projection. Adaptive filtering and blocked, matrix-oriented implementations manage convergence, memory, communication, and computational cost.

  • Filtering role: Chebyshev filtering adaptively approximates the lowest occupied eigenspace of the finite-element Hamiltonian.The filtering acts on a trial subspace using a scaled Hamiltonian whose wanted spectrum lies below −1.
  • ChFSI procedure: ChFSI proceeds through Chebyshev filtering, orthonormalization, and Rayleigh–Ritz projection to obtain the occupied subspace.The SCF iteration then computes a new output electron density.
  • Adaptive filtering: An adaptive filtering strategy repeats ChFSI until the residual norm near the Fermi energy falls below δ = 1×10^-2–5×10^-2.The polynomial degree is typically 20–50 and depends on the Hamiltonian spectrum bound.
  • Implementation: Chebyshev filtering costs O(MN), where M is the discretized Hamiltonian size and N is the number of occupied states.DFT-FE optimizes this kernel using stored cell-level Hamiltonian matrices and BLAS Xgemm operations.
  • Implementation: Filtering wavefunctions in blocks reduces peak memory, with a default block-size range of 300–400.Very small blocks increase memory-access and latency overheads, while very large blocks increase memory usage and reduce batched Xgemm efficiency.

3.4. Cholesky factorization based Gram-Schmidt orthonormalization

DFT-FE uses Cholesky-based Gram-Schmidt orthonormalization with blocked, distributed computations to control memory, communication, and computational costs for large systems. Mixed precision accelerates the dominant overlap and basis-construction steps while preserving SCF accuracy and robustness.

  • Cholesky factorization based Gram-Schmidt orthonormalization: For systems larger than 20,000 electrons, orthonormalization becomes a dominant cubic-scaling cost, motivating the CholGS procedure in DFT-FE.CholGS uses Cholesky factorization based Gram-Schmidt orthonormalization to prevent ill-conditioning of filtered vectors.
  • Parallel implementation aspects: The overlap matrix S is distributed with a 2D cyclic block grid, and only its lower triangular part is computed to reduce memory and work.The blocked computation successively forms submatrices while exploiting Hermiticity.
  • Parallel implementation aspects: Cholesky factorization and inverse-factor computation together contributed about 7% of CholGS wall time for a 61,502-electron system.These operations use ScaLAPACK routines pXpotrf and pXtrtri and were minor relative to the other CholGS steps.
  • Parallel implementation aspects: A two-level blocked construction of the orthonormalized basis reuses wavefunction storage and avoids replicating the full inverse Cholesky factor on every MPI task.Outer blocking acts over locally owned FE nodes, while inner blocking acts over wavefunction columns.
  • Mixed precision approaches: Mixed precision makes the off-diagonal overlap and corresponding basis-construction computations single precision while retaining double precision for diagonal operations.Single-precision collectives also halve the relevant MPI message size in bytes compared with double precision.
  • Mixed precision approaches: Around 2 times faster, mixed precision reduced the dominant CholGS costs for benchmark systems containing 39,900 and 61,502 electrons.The paper also evaluates the resulting SCF convergence accuracy and robustness.

3.5. Rayleigh-Ritz procedure and electron-density computation

The Rayleigh-Ritz procedure is a major large-scale bottleneck, so DFT-FE combines spectrum splitting, partial diagonalization, blocking, and mixed precision to reduce its cost. These strategies retain the relevant fractionally occupied states and reduce projected-Hamiltonian work while maintaining numerical accuracy.

  • Rayleigh-Ritz procedure: Rayleigh-Ritz includes projected-Hamiltonian construction, diagonalization, and subspace rotation, with computational complexities O(MN^2), O(N^3), and O(MN^2), respectively.Together, these operations make Rayleigh-Ritz a significant bottleneck for large-scale problems.
  • Spectrum splitting: Spectrum splitting computes only the Nfr largest eigenstates needed for electron-density evaluation, while fully occupied states are handled separately.For metallic systems, Nfr is typically 10–15% of N and is much smaller for some other systems.
  • Spectrum splitting: Replacing O(MN^2) subspace rotation with O(MNNfr) reduces the dependence on the total number of states when only fractionally occupied states are required.Partial diagonalization is performed for the Nfr largest eigenstates.
  • Mixed precision in RR: Mixed precision computes the fractionally occupied projected-Hamiltonian block in double precision and the other blocks in single precision.Because Nfr is typically less than 15% of N, the double-precision block is a small part of the computation.
  • Mixed precision in RR: Around 2 in computation of the projected Hamiltonian, mixed precision reduced its cost by about a factor of 2 while preserving SCF counts and causing negligible energy and force errors.The benchmark discretization errors were approximately 10^-4 Ha/atom in energy and 10^-4 Ha/Bohr in ionic forces.

3.6. Mixing schemes

DFT-FE treats SCF convergence as a fixed-point problem in electron density and uses Anderson mixing in all reported simulations. More advanced mixing strategies are planned for a future release.

  • Mixing schemes: The SCF iteration is formulated as the fixed-point problem ρ = F[ρ], with electron density mixed to accelerate convergence.F[ρ] represents computation of the occupied eigenspace for a given density.
  • Mixing schemes: Anderson mixing was used in all simulations conducted in the work.The authors state that more advanced mixing strategies are planned for a future release.
  • Mixing schemes: Future releases are intended to implement mixing strategies that provide improved SCF convergence rates independent of system size.This is presented as planned work rather than a capability evaluated in the current simulations.

3.7. Parallelization

DFT-FE combines domain decomposition with wavefunction band parallelization to improve scalability across SCF operations. A measured hybrid strategy outperformed domain decomposition alone at 102,400 MPI tasks.

  • Domain and band parallelization: Domain decomposition partitions the adaptive FE mesh across MPI tasks, with communication restricted primarily to FE nodes on processor boundaries.The localized FE basis has compact support, reducing communication relative to globally coupled representations.
  • Domain and band parallelization: Band parallelization divides wavefunctions into groups distributed across MPI subcommunicators for filtering, orthonormalization, Rayleigh-Ritz, and electron-density computation.Electron-density contributions from band groups require a small accumulation communication.
  • Domain and band parallelization: Increasing band-parallelization groups can significantly raise communication costs because all-to-all exchanges involve large wavefunction, overlap, and projected-Hamiltonian data.The paper therefore combines only a modest amount of band parallelization with domain decomposition.
  • Scalability comparison: At 102,400 MPI tasks, the combined approach P3 achieved 49% efficiency, compared with 41% for domain decomposition alone.P3 used domain decomposition through 51,200 tasks at 73% efficiency, then added two band-parallelization groups.

3.8. Geometry optimization algorithm

DFT-FE performs geometry optimization with a nonlinear conjugate-gradient method and regenerates the adaptive finite-element mesh after each atomic-force update. Future work targets more robust line searches and faster-converging optimization algorithms.

  • DFT-FE uses the Polak–Ribiere-Polyak nonlinear conjugate-gradient algorithm with a secant line search for geometry optimization.
  • Atomic-force relaxation regenerates the finite-element mesh after every geometry update step.
  • Future work plans to replace the secant line search with Brent’s method and implement L-BFGS and FIRE algorithms.The paper notes that L-BFGS and FIRE demonstrate faster convergence than PRP-CG in many cases.

4. Results and discussion

The results establish DFT-FE’s accuracy across pseudopotential and all-electron benchmarks, strong scaling to large MPI task counts, and increasing computational advantages over plane-wave codes as systems grow.

  • Accuracy validation: DFT-FE benchmarks cover periodic and non-periodic pseudopotential and all-electron calculations, with comparisons against QE, ABINIT, exciting, and NWChem.The validation examines energies, ionic forces, and cell stresses across systems ranging from 2,550 to 39,990 electrons for scalability studies.
  • Accuracy validation: DFT-FE agrees excellently with QE for Mg, Mo, and Cu vacancy benchmarks, with differences systematically decreasing as discretization increases.
  • Accuracy validation: All-electron benchmarks show good agreement with exciting and NWChem, including approximately 1 × 10^-4 Ha/atom energy agreement for the Si vacancy system.The Si vacancy comparison also reports approximately 5 × 10^-4 Ha/Bohr ionic-force and 5 × 10^-6 Hartree/Bohr^3 hydrostatic-stress agreement.
  • Parallel scaling: 75% scaling efficiency is achieved for Mg4x4x4 on 4,096 MPI tasks, with a 28-second wall time per SCF iteration.The study attributes this scalability for small-to-medium systems to low communication costs in finite-element Hamiltonian and wavefunction products.
  • Parallel scaling: 49% scaling efficiency is achieved for Mg10x10x10 on 102,400 MPI tasks, with a 237-second wall time per SCF iteration.For Mg8x8x8, combined domain and band parallelization reaches 64,000 MPI tasks at 43% efficiency.
  • Computational efficiency: DFT-FE’s CPU-time advantage begins beyond approximately 3,000 electrons for hcp Mg and 6,000 electrons for bcc Mo.
  • Computational efficiency: Close-to-quadratic scaling extends to 39,990 electrons for Mg, 27,986 for Mo, and 17,537 for Cu nanoparticles.The reported complexity exponents are O(N_e^2.12), O(N_e^2.32), and O(N_e^2.04), respectively.

5. Demonstration of DFT-FE’s capabilities

DFT-FE is demonstrated on ionic relaxation, molecular dynamics, band structures, and defect states. These examples show validation against established codes, stringent energy conservation, and access to spin-polarized defect physics.

  • Geometry optimization: DFT-FE performs ionic relaxation of tris(bipyridine) ruthenium and validates energies and forces against QE.The 61-atom complex is evaluated with energy and force discretization errors converged to approximately 10^-4 Ha per atom and 10^-4 Ha/Bohr.
  • Implementation diagnostics: DFT-FE’s wall-time breakdown separates electrostatic solves, Chebyshev filtering, orthogonalization, Rayleigh-Ritz, density computation, and other operations.
  • Molecular dynamics: An NVE simulation of 108-atom fcc Al uses velocity Verlet integration with a 0.5 fs timestep and initial ionic temperature of 1,500 K.
  • Molecular dynamics: The Al molecular-dynamics run has total-energy standard deviation 1.482 × 10^-5 Ha/atom and energy drift 1.8045 × 10^-9 Ha/atom-fs.The reported mean total energy is −2.0948 Ha/atom over 800 fs.
  • Band-structure calculations: DFT-FE computes silicon band structures by extracting Kohn-Sham eigenvalues along high-symmetry Brillouin-zone lines using an 8 × 8 × 8 Monkhorst-Pack grid.The calculation uses an ONCV pseudopotential and PBE functionals, with comparison to QE.
  • Defect-state calculations: A 4 × 4 × 4 silicon supercell with a vacancy produces spin-polarized midgap states, including majority-spin energy 0.54700 eV and minority-spin energy 0.79010 eV.The top of the valence band is used as the energy reference.

6. Conclusions

DFT-FE combines adaptive higher-order finite elements with locality-aware algorithms and careful implementation to enable accurate, scalable large-scale DFT. Validation shows excellent agreement with established codes and substantial CPU-time and wall-time advantages over QE.

  • DFT-FE supports pseudopotential and all-electron calculations on periodic, non-periodic, and semi-periodic systems, with adaptive spectral finite-element bases up to 12th order.Spatial adaptivity can use user-defined mesh parameters or automatic refinement based on local error estimates.
  • Chebyshev filtering, Cholesky-based orthonormalization, Rayleigh-Ritz projection, mixed precision, and spectrum splitting reduce computational prefactors, communication costs, and latency.Spectrum splitting reduces the prefactors of CholGS and Rayleigh-Ritz by around 2 and 3, respectively, for large systems.
  • DFT-FE agrees excellently with QE, exciting, and NWChem for benchmark ground-state energies, ionic forces, and stresses.For pseudopotential calculations, differences from QE were significantly lower than the codes’ discretization errors.
  • 4.5−12× lower CPU-times were obtained than QE for periodic systems containing 10,659−20,470 electrons.DFT-FE was more efficient than QE beyond 3,000 electrons.
  • DFT-FE provides a practical capability for large-scale pseudopotential calculations reaching 100,000 electrons and can also perform all-electron calculations in the same framework.The framework is intended for generic material systems with arbitrary boundary conditions and complex geometries.
Loading 1903.10959v3…