SPARC: Accurate and efficient finite-difference formulation and parallel implementation of Density Functional Theory: Extended systems

1. Introduction

• Plane-wave discretization is popular for its complete orthonormal basis, spectral convergence, efficient FFT convolutions, and effective preconditioning, but periodicity and non-locality constrain its use for localized systems and large-scale computing.The text also notes that these limitations make linearly scaling methods impractical.Source: §1 · p. 1
• Real-space approaches include finite differences, finite elements, wavelets, periodic sinc functions, splines, and maximum-entropy basis functions, yet plane-wave implementations remain the practical method of choice.The text attributes this to their superior chemical-accuracy performance on commonly available computational resources.Source: §1 · pp. 1–2
• Finite differences yield a standard eigenvalue problem with relatively small spectral width and readily support high-order approximations, but can produce non-monotonic energy and force convergence and egg-box effects.The egg-box effect can affect structural relaxations and molecular dynamics.Source: §1 · p. 2
• The earlier SPARC component addressed isolated systems using local electrostatics, Chebyshev-filtered self-consistent field iteration, and a reformulated non-local force component.For isolated systems, SPARC was reported to outperform ABINIT, PARSEC, and OCTOPUS in achieving desired electronic-structure accuracy.Source: §1 · p. 2
• SPARC shows high energy and force convergence to plane-wave results, negligible egg-box effects and molecular-dynamics drift, and plane-wave-like scaling with a smaller prefactor for systems up to thousands of electrons.The work also reports accurate ground-state properties and exponential convergence with vacuum size for slabs and wires.Source: §1 · p. 2
• Extended systems require careful electrostatic reformulation, Bloch-periodic orbitals, and Brillouin-zone integration, motivating a separate investigation of DFT accuracy and efficiency.The paper targets crystals, surfaces, and wires.Source: §1 · p. 2

2. Density Functional Theory (DFT)

• The electron density is expressed using orbitals and their occupations, with the Bloch wavevector and Brillouin-zone volume average entering the extended-system description.The orbitals and occupations are denoted Ψ and g, respectively.Source: §2 · p. 3
• The unit-cell free energy combines kinetic, exchange-correlation, non-local pseudopotential, electrostatic, and entropy contributions.The entropy term arises from partial orbital occupations.Source: §2 · pp. 3–4
• The exchange-correlation energy has no known exact form, so approximations such as the Local Density Approximation are used.The LDA exchange-correlation energy per particle combines exchange and correlation contributions for a uniform electron gas.Source: §2 · p. 3
• The non-local pseudopotential energy uses the Kleinman-Bylander representation, with coefficients and projection functions defined using isolated-atom pseudowavefunctions and pseudopotentials.The angular-momentum indices l and m denote azimuthal and magnetic quantum numbers.Source: §2 · p. 3
• Electrostatic energy comprises electron-density Hartree energy, electron-nucleus interaction, and nucleus-nucleus repulsion, including periodic images.The entropy contribution originates from partial orbital occupations.Source: §2 · p. 4
• The ground state is obtained through a variational problem, and its electronic state and Hellmann-Feynman forces are recomputed for each nuclear configuration during geometry optimization or molecular dynamics.The system is specified by a unit cell containing N atoms and Ne valence electrons.Source: §2 · pp. 3–4

3. Formulation and implementation

• SPARC replaces the O(N^2) direct real-space evaluation of electrostatic energy, whose components diverge in extended systems, with a local formulation.The formulation uses the electrostatic potential and total nuclear pseudocharge density.Source: §3 · p. 4

3.1. Finite-difference discretization

• Central finite differences approximate the Laplacian and gradient on a cuboid grid, with trapezoidal spatial integration and periodic or zero-Dirichlet boundary conditions.SPARC uses PETSc data structures and routines for its implementation.Source: §3 · p. 6, §3.1 · pp. 7–8

3.2. Pseudocharge density generation and self energy calculation

• The finite-difference Laplacian generates pseudocharge densities, whose exponential decay away from each nucleus permits truncation.The truncation region is chosen to satisfy a prescribed charge tolerance.Source: §3.2 · pp. 8–9

3.3. Brillouin zone integration

• Brillouin-zone integration approximates volume-averaged integrals using weighted sampling nodes, with the node choice defining the sampling scheme.Source: §3.3 · pp. 9–10

3.4. Electron density calculation

• CheFSI computes approximations to eigenvectors for the lowest Ns eigenvalues at each k-point, using Chebyshev filtering to remove unwanted components.The authors motivate this eigensolver by its low orthogonalization and memory costs relative to alternatives; the first SCF iteration performs multiple CheFSI steps, while later iterations perform one.Source: §3.4 · pp. 10–11

3.5. Free energy calculation

• The discrete ground-state free-energy evaluation scales as O(Nd), and therefore O(N) with respect to atom count.It uses discretized spatial and Brillouin-zone integrals and includes self-energy and overlapping-pseudocharge corrections.Source: §3.5 · p. 11

3.6. Atomic forces calculation

• The discrete atomic force separates into local, overlapping-pseudocharge correction, and non-local components; transferring projector derivatives to orbitals improves force accuracy at a given discretization.This introduces an approximation because the finite-difference gradient is approximate, but the authors report that smoother orbitals offset it.Source: §3 · p. 6, §3.6 · p. 11
• The local force component is calculated using Eqn. 39 and outlined in Algorithm 2.Source: §3.6 · p. 12
• Approximating the spatial integral with Eqn. 24 yields the discrete non-local force component.Source: §3.6 · p. 12
• Algorithm 3 outlines calculation of the non-local force component.Source: §3.6 · pp. 12–13
• For each atom, the non-local force calculation uses a cuboid centered on that atom, chosen over a sphere for simplicity and efficiency in the Euclidean finite-difference discretization.The cuboid has side lengths 2rc.Source: §3.6 · p. 12
• Within the processor overlap region, radial projector components are interpolated onto the finite-difference grid using cubic splines.Source: §3.6 · p. 12
• Processor contributions to YJnlm and WJnlm are summed, and the resulting quantities are used to calculate the non-local force with Eqn. 40.Source: §3.6 · p. 12
• Atomic-force calculation scales as O(N) with the number of atoms.Source: §3.6 · p. 13

4.1. Convergence with discretization

• Mesh refinement produced systematic convergence to plane-wave results, with average rates of O(h^10) for energy and O(h^9) for forces.The study used lithium hydride, silicon, and gold; chemical accuracy was readily attained without double-grid or high-order integration measures.Source: §4.1 · p. 14, Figure 2 · p. 14

4.2. Convergence with vacuum size for slabs and wires

• For slabs and wires, SPARC’s energies and forces converged exponentially with vacuum size; around 7 Bohr sufficed for chemical accuracy, while ABINIT’s error stagnated.The ABINIT stagnation was attributed to spurious interactions between periodic images.Source: §4.2 · p. 15, Figure 3 · p. 15

4.3. Bulk properties

• For bulk silicon, SPARC and ABINIT energy curves were practically indistinguishable; equilibrium lattice constant, energy, and bulk modulus agreed within 0.003 Bohr, 1×10^-5 Ha/atom, and 0.18 GPa.Their band structures were nearly identical, with HOMO, LUMO, and bandgap agreeing within 7×10^-5 Ha, 1×10^-5 Ha, and 6×10^-5 Ha, respectively.Source: §4.3 · pp. 15–16, Figure 4 · p. 16, Figure 5 · p. 16

4.4. Slab and wire properties

• For the aluminum slab, surface energies agreed within 0.001 J/m2 and maximum force differences were within 3×10^-4 Ha/Bohr; silicon-wire DOS curves were practically indistinguishable.For the wire, energy and force differences were within 7×10^-5 Ha/atom and 3×10^-4 Ha/Bohr, respectively.Source: §4.4 · p. 17, Figure 6 · p. 17, Figure 7 · p. 18

4.5. Geometry optimization

• Computed forces matched the energy derivative without noticeable egg-box effects at chemical-accuracy meshes; vacancy formation energies agreed within 7×10^-4 Ha and improved to 9×10^-6 Ha on mesh refinement.For the vacancy systems, relaxed atomic positions differed by at most 1.6×10^-3 Bohr initially and 1.3×10^-4 Bohr on the refined meshes.Source: §4.5 · pp. 18–19, Figure 8 · p. 19, Figure 9 · p. 20

4.6. Molecular dynamics

• In the NVE simulation, the linear-fit total-energy drift was 1 × 10^-4 Ha/atom-ps, and no significant drift was observed overall.9The mean total energy was −2.07842 Ha/atom, with a standard deviation of 1.6×10^-4 Ha/atom.Source: §4.6 · p. 20, Figure 10 · p. 21

4.7. Scaling and Performance

• SPARC showed similar strong-scaling trends to ABINIT, with prefactors smaller by up to 4×; weak-scaling speedups were 2.2–6×, and minimum wall times were over 2.5× lower, approximately 5× for the largest system.The weak-scaling tests covered 321–4116 electrons.Source: §4.7 · p. 22

5. Concluding Remarks

• SPARC combines a finite-difference representation, Chebyshev-filtered self-consistent field iteration, and reformulated electrostatics and non-local forces to evaluate energies and atomic forces.The framework targets extended systems including crystals, slabs, and wires.Source: §5 · p. 23
• SPARC shows high mesh-size convergence toward plane-wave results, exponential vacuum-size convergence for slabs and wires, negligible egg-box effects and molecular-dynamics drift, and similar scaling to optimized plane-wave codes with a smaller prefactor.The scaling comparison covers systems with up to thousands of electrons; SPARC also significantly outperforms other finite-difference DFT packages.Source: §5 · p. 23
• Band parallelization, scalable parallel diagonalization of the subspace eigenproblem, further code optimization, and removal of PETSc and MKL dependencies are identified as future improvements.The authors expect the first two changes to improve strong scaling and enable significantly larger systems, and removal of external-package dependencies to improve portability and usability.Source: §5 · p. 23

A. Electrostatic correction for overlapping pseudocharge densities

• Overlapping pseudocharge densities create a disagreement with the conventional repulsive-energy treatment, which treats nuclei as point charges.The paper states that its electrostatic formulation does not make this distinction and introduces a correction to reestablish agreement.Source: A. Electrostatic correction for… · p. 28
• The reference pseudocharge density and each nucleus's spherically symmetric, compactly supported reference charge density are used to define the reference potential.The latter generates the potential associated with that nucleus.Source: A. Electrostatic correction for… · p. 28
• The repulsive-energy correction is discretized using integrals approximated by the trapezoidal rule.The discrete correction is given as the form of Eqn. 46.Source: A. Electrostatic correction for… · p. 28
• The correction to atomic forces is also specified, with the summation over an atom and its periodic images.The reference potential is the one previously used to generate neutralizing densities in all-electron calculations.Source: A. Electrostatic correction for… · pp. 28–29

B. Pseudopotential parameters

• The l = 0 component is chosen to be local in all calculations.The listed parameters concern the different angular momentum components of the Troullier-Martins pseudopotential.Source: B. Pseudopotential parameters · p. 29
• The listed pseudopotential parameters distinguish angular momentum components, with the l = 0 component used as local throughout the calculations.The table is identified as covering non-local projector cutoff radii.Source: B. Pseudopotential parameters · p. 29, Table 2 · p. 29
• Table 2 lists cutoff radii for non-local projectors within the Troullier-Martins pseudopotential.The supplied passage identifies the table's scope but does not provide its values.Source: Table 2 · p. 29

C. Discrete pseudocharge density properties

• The discrete pseudocharge density has infinite extent because the Laplacian is approximated by finite differences, despite compact support of the continuous density.The continuous density is supported within the local-pseudopotential cutoff radius around its atom.Source: C. Discrete pseudocharge density… · p. 29
• At h = 0.5 Bohr, the normalized net-enclosed-charge error decays exponentially with pseudocharge radius, allowing truncation at a finite radius without significant accuracy loss.Figure 12 plots the error as a function of pseudocharge radius.Source: C. Discrete pseudocharge density… · p. 29, Figure 12 · p. 30
• The radius needed for εb = 10^-8 decreases as mesh size is refined and tends to zero as h approaches zero.The slight non-monotonicity occurs because SPARC chooses the radius as a multiple of h.Source: C. Discrete pseudocharge density… · p. 30, Figure 13 · p. 30

• For 3 × 3 × 3 FCC aluminum cells with a vacancy, SPARC has the best energy and force convergence among the finite-difference codes; PARSEC and OCTOPUS energy errors use an h = 0.201 Bohr baseline.The comparison plots convergence against mesh size; the PARSEC and OCTOPUS energy baseline differs from ABINIT.Source: D. Comparison of SPARC with other… · pp. 30–31, Figure 14 · p. 31
• In strong scaling, SPARC’s minimum wall time is lower by factors of 12.2, 93, and 5.8 than PARSEC, OCTOPUS, and ABINIT, respectively.The authors attribute SPARC’s advantage over PARSEC to higher strong-scaling efficiency, and over OCTOPUS to a smaller prefactor and superior strong scaling.Source: D. Comparison of SPARC with other… · p. 31
• SPARC’s local electrostatic reformulation outperforms conventional Fourier-based approaches, with its electrostatics time 2.7× lower than OCTOPUS and the speedup reaching 62.6× on 384 cores.Source: D. Comparison of SPARC with other… · p. 31
• SPARC achieves minimum wall-time speedups of up to 12.2×, 72.3×, and 5.55× versus PARSEC, OCTOPUS, and ABINIT, respectively.The comparison excludes the first SCF iteration and covers aluminum vacancy-cell systems; Table 3 presents the minimum wall times.Source: D. Comparison of SPARC with other… · p. 32, Table 3 · p. 32