PaperlayerView exact revision on arXiv

Source-linked AI summary

SPARC: Accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory: Isolated clusters

Swarnava Ghosh, Phanish Suryanarayana

arXiv:1603.04334v3physics.comp-ph

TL;DR

DFT’s computational cost and plane-wave methods’ limitations for isolated clusters motivate a practical real-space alternative. SPARC uses finite differences with reformulated electrostatics and forces and filtered self-consistent-field iteration; the paper reports accurate energies and forces and parallel scaling similar to optimized plane-wave codes, with a smaller prefactor.

  • Problem

    DFT remains computationally demanding, while plane-wave methods require artificial supercells for clusters and real-space implementations have struggled to match their accuracy and performance on modest resources.

    Source: §1 · pp. 1–2

  • Method

    SPARC formulates DFT for isolated clusters using finite differences, local electrostatics, Chebyshev-filtered self-consistent-field iteration, and a reformulation of the non-local atomic force.

    Source: §1 · p. 3, §3 · p. 7, §3.2 · p. 9

  • Results

    SPARC reports exponential domain-size convergence, mesh-refined convergence toward plane-wave results, consistent forces without noticeable egg-box effects, accurate ground-state properties, and similar parallel scaling with a significantly smaller prefactor.

    Source: §1 · p. 3, §4.1 · p. 16, §4.4 · pp. 20–21

  • Takeaways & Limitations

    SPARC is presented as an attractive alternative to plane-wave codes for practical DFT simulations of isolated clusters.

    Source: §1 · p. 3

    The examples are restricted to about 1000 atoms, and the serial subspace eigenvalue solve is expected to dominate for systems with tens of thousands of electrons.

    Source: §5 · p. 22

Abstract

from arXiv · show

As the first component of SPARC (Simulation Package for Ab-initio Real-space Calculations), we present an accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory (DFT) for isolated clusters. Specifically, utilizing a local reformulation of the electrostatics, the Chebyshev polynomial filtered self-consistent field iteration, and a reformulation of the non-local component of the force, we develop a framework using the finite-difference representation that enables the efficient evaluation of energies and atomic forces to within the desired accuracies in DFT. Through selected examples consisting of a variety of elements, we demonstrate that SPARC obtains exponential convergence in energy and forces with domain size; systematic convergence in the energy and forces with mesh-size to reference plane-wave result at comparably high rates; forces that are consistent with the energy, both free from any noticeable `egg-box' effect; and accurate ground-state properties including equilibrium geometries and vibrational spectra. In addition, for systems consisting up to thousands of electrons, SPARC displays weak and strong parallel scaling behavior that is similar to well-established and optimized plane-wave implementations, but with a significantly reduced prefactor. Overall, SPARC represents an attractive alternative to plane-wave codes for practical DFT simulations of isolated clusters.

1. Introduction

DFT offers high accuracy relative to cost but remains computationally demanding, while plane-wave methods are poorly suited to isolated clusters and real-space methods have struggled to match their practical performance. SPARC presents a finite-difference approach for isolated clusters that targets accurate energies and forces alongside competitive parallel scaling.

1. Introduction

  • DFT’s accuracy-to-cost ratio has driven broad use, but orthogonality and large basis sets produce cubic scaling and substantial global communication.These costs restrict accessible system sizes, particularly for ab initio molecular dynamics requiring tens to hundreds of thousands of Kohn–Sham solves.Source: §1 · p. 1
  • Plane-wave methods require artificial supercell periodicity for clusters and defects, and their non-locality complicates linear scaling and distributed-memory parallelization.The spherical spectral scheme for isolated clusters inherits these characteristics.Source: §1 · p. 2
  • Real-space methods offer flexible boundary conditions and routes to linear scaling and scalable computing, yet have not consistently matched optimized plane-wave accuracy and performance on modest resources.Plane-wave codes also provide substantially more functionality after a longer development history.Source: §1 · p. 2
  • Finite differences yield a standard eigenproblem, relatively small spectral width, compact Laplacian representation, and straightforward high-order approximations.Their lack of an underlying basis and associated variational structure can cause non-monotonic energy and force convergence; lower-order integration can also produce egg-box effects.Source: §1 · pp. 2–3
  • SPARC combines local electrostatics, Chebyshev-filtered self-consistent-field iteration, and a reformulated non-local force to compute accurate energies and atomic forces for isolated clusters.The paper reports exponential convergence with domain size, convergence toward plane-wave results with mesh refinement, forces consistent with energy and without noticeable egg-box effects, accurate ground-state properties, and similar scaling with a smaller prefactor.Source: §1 · p. 3

2. Density Functional Theory (DFT)

This section formulates finite-temperature DFT for isolated systems as a variational problem over orbitals, occupations, and nuclear positions. It describes the energy contributions, electron density, and approximations used in the formulation.

2. Density Functional Theory (DFT)

  • The DFT free energy combines kinetic, exchange-correlation, non-local pseudopotential, electrostatic, and electronic-entropy contributions.The ground state is posed as a variational problem, and geometry optimization requires computing the electronic ground state for each nuclear configuration.Source: §2 · pp. 3–5
  • The electron density depends on the orbitals and their occupations, while the free-energy formulation includes nuclear positions and electronic temperature.The system is described as isolated, with N atoms and Ne valence electrons, neglecting spin.Source: §2 · p. 3
  • The work employs the Local Density Approximation for exchange-correlation energy because its exact form is unknown.The section also identifies the Generalized Gradient Approximation as another popular approximation.Source: §2 · p. 4
  • The non-local pseudopotential uses the separable Kleinman–Bylander form, with projection functions defined using isolated-atom pseudowavefunctions and angular-momentum-dependent pseudopotentials.The local pseudopotential component is typically chosen as one angular-momentum-dependent component.Source: §2 · p. 4
  • Electrostatic energy is decomposed into electron-density Hartree energy, electron–nucleus interaction, and nucleus–nucleus repulsion.Electronic entropy accounts for partial orbital occupations using a dependence appropriate for Fermions.Source: §2 · p. 4

3. Formulation and implementation

This section presents SPARC’s real-space finite-difference formulation and parallel implementation for DFT simulations of isolated clusters, covering localized electrostatics, ground-state and electron-density calculations, and accurate energy and force evaluation. It also describes interpolation and parallel-reduction procedures and reports linear scaling of force evaluation with atom count.

3. Formulation and implementation

  • Local electrostatics replaces a non-local energy evaluation that scales as O(N^2) in atom count and requires substantial interprocessor communication.The electronic ground state is found with SCF iterations, while the orthonormal Kohn–Sham orbital calculation scales asymptotically as O(N^3) in atom count.Source: §3 · pp. 5–6
  • SPARC transfers the position derivative of non-local projectors to spatial derivatives of the smoother orbitals to improve atomic-force accuracy.The force expression also includes local and pseudocharge-overlap electrostatic correction terms.Source: §3 · p. 7
  • SPARC uses PETSc data structures, stores the discrete Laplacian in compressed-row format, applies non-local pseudopotentials matrix-free, and stores orbitals as dense-matrix columns.The local pseudocharge-generation and self-energy computations scale as O(N) in atom count and have good weak and strong parallel scalability, although atom assignment makes the overall procedure formally slightly worse than O(N).Source: §3 · p. 7, §3.1 · p. 8, §3.2 · pp. 9–10

3.1. Finite-difference discretization

  • Finite-difference Laplacian and gradient approximations have O(h^2no) accuracy, while midpoint integration makes the discrete free energy consistent with the discrete Kohn–Sham equations.The scheme uses zero Dirichlet values outside the grid and a uniform finite-difference mesh.Source: §3.1 · p. 8

3.2. Pseudocharge density generation and self energy calculation

  • Discrete pseudocharge densities have exponentially decaying tails that permit truncation, and overlap corrections rapidly converge to zero with mesh refinement when continuous pseudocharges do not overlap.Higher-order finite differences produce smaller truncation radii even when pseudopotentials are not smooth at their cutoff.Source: §3.2 · p. 9

3.3. Electrostatic potential calculation

  • A discrete multipole charge correction imposes electrostatic boundary conditions because the potential decays algebraically, making zero Dirichlet boundaries alone inaccurate, especially for charged or dipolar systems.With multigrid, the potential calculation is stated to scale as O(N) in atom count.Source: §3.3 · pp. 10–11

3.4. Electron density calculation

  • CheFSI computes approximations to the lowest eigenpairs by filtering orbitals with Chebyshev polynomials, then projecting, solving a smaller eigenproblem, and rotating the subspace.Filtering damps eigencomponents above a chosen cutoff; the method is selected for lower orthogonalization and memory costs. In later SCF iterations, the rotated orbitals provide the initial guess.Source: §3.4 · pp. 11–13
  • The electron-density subsection reports O(N^3) scaling with respect to the number of atoms.Source: §3.4 · p. 13

3.5. Free energy calculation

  • The discrete free-energy expression is obtained by approximating the integrals in Eqn. 16 with the integration rule in Eqn. 22.Source: §3.5 · p. 13
  • Free-energy evaluation scales as O(Nd), hence O(N) with atom count, and is performed at every SCF step despite being needed only after the ground state is determined.Source: §3.5 · p. 13

3.6. Atomic forces calculation

  • The discrete atomic force comprises electrostatic-correction, local, and non-local components.Source: §3.6 · p. 13
  • The local force procedure interpolates the potential onto the finite-difference grid with cubic splines in the overlap region and sums processor contributions for all atoms.Source: §3.6 · pp. 13–14
  • For non-local forces, SPARC uses a cube of side 2r_c for each atom's cutoff region, interpolates with cubic splines, and sums processor contributions.The cube is chosen over a sphere for simplicity and efficiency in the Euclidean finite-difference discretization.Source: §3.6 · p. 14
  • Atomic-force evaluation scales as O(N) with respect to the number of atoms.Source: §3.6 · p. 15

4. Examples and Results

SPARC’s isolated-cluster calculations were tested for domain- and mesh-size convergence, ground-state accuracy, and parallel performance against ABINIT. The results show accurate energies and forces alongside competitive scaling with lower wall-time prefactors.

4.1. Convergence with domain size

  • At a 12 × 12 × 12 Bohr domain, CO and H2O already reached chemical accuracy in energy and forces, which converged exponentially as the domain grew.The tests used deliberately asymmetric molecular placement to exaggerate finite-domain effects.Source: §4.1 · p. 16

4.2. Convergence with spatial discretization

  • Mesh refinement produced systematic convergence toward ABINIT, with average rates of approximately O(h^9) for energy and O(h^8) for forces.The comparisons covered Pt13, Au13, and C13H33NO2Si3; energy and force convergence rates were comparable without double-grid or high-order integration techniques.Source: §4.2 · p. 17, Figure 4 · p. 18

4.3. Ground state properties

  • For H2, N2, and O2, computed forces agreed with energy derivatives without a noticeable egg-box effect, and vibrational frequencies differed from ABINIT by at most 8 cm−1.Geometry optimization also agreed with ABINIT to within 0.0007 Ha/atom in energy and 0.001 Bohr in equilibrium bond length.Source: §4.3 · pp. 18–20, Table 1 · p. 19, Table 2 · p. 20

4.4. Scaling and performance

  • SPARC and ABINIT showed similar strong-scaling trends for Si275H172, with SPARC outperforming ABINIT by up to 7×.No further wall-time reduction was observed after approximately 700 cores for SPARC and 600 cores for ABINIT.Source: §4.4 · p. 20, Figure 7 · p. 21
  • SPARC’s lower performance prefactor was attributed primarily to more efficient Hamiltonian application in the finite-difference discretization.For larger systems, ABINIT’s nonlocal operations became dominant; excluding those operations, SPARC still achieved speedups of 4.2× for Si29H36 and 1.8× for Si275H172 in weak scaling.Source: §4.4 · p. 21
  • SPARC achieved wall times more than 6.8× lower than ABINIT across the reported systems, with similar strong- and weak-scaling trends and weak-scaling speedups of 3–7×.For Si525H276, SPARC required approximately 7.0× less wall time; the weak-scaling comparison covered clusters from 152 to 3768 electrons.Source: §4.4 · pp. 20–22

5. Concluding Remarks

SPARC combines finite differences with reformulated electrostatics and atomic forces to calculate energies and forces for isolated clusters. The paper reports convergence and ground-state accuracy alongside plane-wave-like parallel scaling at a smaller prefactor.

5. Concluding Remarks

  • SPARC achieves exponential energy and force convergence with domain size, mesh-size convergence comparable to reference plane-wave results, and forces consistent with energy without noticeable egg-box effects.The reported accurate ground-state properties include equilibrium energies, geometries, and vibrational spectra across examples with light and heavy elements.Source: §5 · p. 22
  • For systems up to thousands of electrons, SPARC’s weak and strong parallel scaling is similar to optimized plane-wave codes but has a significantly smaller prefactor.Source: §5 · p. 22
  • The examples were restricted to approximately 1000 atoms for thorough accuracy and efficiency analysis, although the authors also report a Si3145H876 simulation.Source: §5 · p. 22
  • The serial subspace eigenvalue problem is expected to become the dominant cost for systems with tens of thousands of electrons; extending SPARC to periodic systems is future work.Source: §5 · p. 22

A. Electrostatic correction for overlapping pseudocharge densities

The electrostatic formulation treats nuclei as point charges for repulsion, unlike its treatment of pseudocharge densities, requiring corrections when those densities overlap. The authors use corrections chosen for simplicity, accuracy, and efficiency in scalable computing.

A. Electrostatic correction for overlapping pseudocharge densities

  • The electrostatic formulation treats nuclei as point charges for repulsion, so overlapping pseudocharge densities produce disagreement with convention unless corrected.Source: A. Electrostatic correction for… · p. 28
  • SPARC uses the corrections in Eqns. 48 and 50 instead of overlap-only corrections that require neighbor lists to be updated at every relaxation step.The authors cite simplicity, accuracy, and efficiency for scalable high-performance computing as reasons for this choice.Source: A. Electrostatic correction for… · p. 29
  • At distances around 3 Bohr, force corrections remain approximately 0.01 Ha/Bohr, supporting their inclusion in structural relaxation and molecular dynamics.Figure 8 plots repulsive-energy and corresponding atomic-force corrections against interatomic distance for N2 and O2.Source: A. Electrostatic correction for… · p. 29, Figure 8 · p. 29

B. Pseudopotential parameters

This subsection identifies the pseudopotential cutoff radii and local-component choice used in the simulations.

B. Pseudopotential parameters

C. Properties of the discrete pseudocharge density

The continuous pseudocharge density has compact support, whereas its discrete form has infinite extent but decays exponentially. A finite truncation radius can be selected to retain accuracy, and the required radius decreases as the mesh is refined.

C. Properties of the discrete pseudocharge density

D. Comparison of SPARC with other finite-difference DFT codes

SPARC is compared with finite-difference and plane-wave DFT codes on mesh convergence and parallel performance. It achieves high convergence rates and lower reported minimum wall times, with the timing comparison excluding the first SCF iteration.

D. Comparison of SPARC with other finite-difference DFT codes

  • Energy converges as O(h^8.85), O(h^7.15), and O(h^8.64) for SPARC, PARSEC, and OCTOPUS; force rates are O(h^9.62), O(h^10.05), and O(h^10.09).The rates are comparable, while PARSEC has a noticeably larger prefactor, particularly for atomic forces.Source: D. Comparison of SPARC with other… · p. 31, Figure 11 · p. 32
  • In weak scaling, SPARC is 3.9× and 5.76× faster than OCTOPUS on Si275H172 and Si525H276, respectively, with the authors attributing this to more RMM-DIIS iterations.OCTOPUS spends more than 70% of its time in GRAM_SCHMIDT, with the share increasing with system size; the authors suggest CheFSI is superior to its default conjugate-gradients eigensolver for real-space DFT.Source: D. Comparison of SPARC with other… · p. 32
  • Strong-scaling minimum wall time is lower by factors of 21 (221 at h=0.3 Bohr in PARSEC), 14, and 6.8 versus PARSEC, OCTOPUS, and ABINIT, respectively.The strong-scaling system is Si275H172, and the first SCF iteration is excluded.Source: D. Comparison of SPARC with other… · p. 32, Figure 12 · p. 33
  • SPARC achieves up to 21× speedup over PARSEC (164× at h=0.3 Bohr), 15× over OCTOPUS, and 7× over ABINIT in minimum wall time.The comparison excludes the first SCF iteration.Source: Table 5 · p. 33, D. Comparison of SPARC with other… · p. 33
Loading 1603.04334v3…