PaperlayerView exact revision on arXiv

Source-linked AI summary

Accurate and efficient linear scaling DFT calculations with universal applicability

Stephan Mohr, Laura E. Ratcliff, Luigi Genovese, Damien Caliste, Paul Boulanger, Stefan Goedecker, Thierry Deutsch

arXiv:1501.05884v1cond-mat.mtrl-sci

TL;DR

Standard Kohn–Sham DFT scales cubically, while existing linear-scaling approaches may limit accuracy or require extensive parameter tuning. BigDFT uses localized support functions with a Daubechies wavelet basis; the authors report accurate, broadly applicable linear scaling with standard settings and moderate computational resources.

  • Problem

    Standard Kohn–Sham DFT scales cubically with system size, and linear-scaling implementations may limit accuracy or require tedious parameter tuning.

    Source: §I · pp. 1–2

  • Method

    BigDFT represents the density matrix using support functions optimized in situ within a Daubechies wavelet basis.

    Source: §I · pp. 1–2

  • Results

    The approach is reported to provide accurate, universally applicable linear-scaling DFT calculations with standard input parameters and moderate computational resources.

    Source: §I · p. 2, §VI · p. 9

  • Takeaways & Limitations

    The authors say standard settings avoid tedious fine-tuning and simplify use of the code by end-users.

    Source: §VI · p. 9

Abstract

from arXiv · show

Density Functional Theory calculations traditionally suffer from an inherent cubic scaling with respect to the size of the system, making big calculations extremely expensive. This cubic scaling can be avoided by the use of so-called linear scaling algorithms, which have been developed during the last few decades. In this way it becomes possible to perform ab-initio calculations for several tens of thousands of atoms or even more within a reasonable time frame. However, even though the use of linear scaling algorithms is physically well justified, their implementation often introduces some small errors. Consequently most implementations offering such a linear complexity either yield only a limited accuracy or, if one wants to go beyond this restriction, require a tedious fine tuning of many parameters. In our linear scaling approach within the BigDFT package, we were able to overcome this restriction. Using an ansatz based on localized support functions expressed in an underlying Daubechies wavelet basis -- which offers ideal properties for accurate linear scaling calculations -- we obtain an amazingly high accuracy and a universal applicability while still keeping the possibility of simulating large systems with only a moderate demand of computing resources.

I. INTRODUCTION

Conventional Kohn–Sham DFT scales cubically, limiting system size, while linear scaling relies on localized quantities. BigDFT combines in-situ optimized support functions with Daubechies wavelets to target accurate, broadly applicable calculations with standard settings and moderate resources.

I. INTRODUCTION

  • Cubic scaling from orthonormalizing extended Kohn–Sham orbitals limits standard DFT to a few hundred or, at most, a few thousand atoms.The cost arises from orbital orthonormality and matrix diagonalization.Source: §I · p. 1
  • Density-matrix elements decay exponentially with distance in insulators and metals at finite temperature, so thresholding small elements can yield linear scaling.The density matrix replaces extended orbitals with a quantity whose nonzero elements can grow linearly with system size.Source: §I · p. 1
  • Daubechies wavelets provide an orthogonal, compact-support basis with reasonable precision at moderate grid spacing; their multiresolution properties permit finer meshes near atoms.These properties make wavelets suitable for representing localized functions while keeping the number of basis functions relatively small.Source: §I · pp. 1–2
  • BigDFT optimizes a smaller set of support functions in situ within predefined localization regions, adapting them to their chemical environment.This approach reduces the subspace used to represent the density matrix.Source: §I · pp. 1–2
  • The BigDFT approach targets accurate, universally applicable linear-scaling DFT without tedious parameter or basis fine-tuning, using default values for reasonable results.The authors say the approach requires only a small amount of computational resources.Source: §I · p. 2
  • About 10 meV/atom absolute energy accuracy is reported for a large variety of isolated systems, while in-situ support functions enable calculations on thousands of atoms with moderate memory and CPU time.The paper proposes CPU minutes per atom as a way to evaluate computational resources for a full DFT calculation.Source: §I · p. 2

II. THEORY AND ALGORITHM

The algorithm represents the density matrix with support functions and a density kernel, using localization and sparsity to achieve linear scaling. It alternates between optimizing support functions and the kernel until overall convergence.

A. Overview of the algorithm

  • The band-structure energy is the central Kohn–Sham quantity calculated using the density-matrix formalism.The Hamiltonian is defined from the Kohn–Sham potential, which includes external, Hartree, and exchange-correlation terms; BigDFT uses norm-conserving pseudopotentials.Source: A. Overview of the algorithm · p. 2
  • A separable density-matrix representation uses support functions and a density kernel, with the charge density given by the density matrix's diagonal.The support functions are expanded in an underlying wavelet basis.Source: A. Overview of the algorithm · pp. 2–3
  • Localization cutoffs for support functions and the density kernel make the support functions localized and the kernel sparse, conditions needed for linear scaling.Kernel elements are set to zero when support-function centers are farther apart than the sum of their kernel cutoffs.Source: A. Overview of the algorithm · p. 2
  • Support functions are optimized with the kernel fixed, then the kernel is optimized self-consistently with fixed support functions; the two inner loops repeat until overall convergence.The first inner loop uses a fixed potential and Hamiltonian, so the algorithm advises exiting it even before its convergence threshold is reached.Source: A. Overview of the algorithm · p. 3, Figure 1 · p. 3

B. Optimization of the support functions and density kernel

The approach optimizes support functions using a target that balances accurate band-structure energy with localization, then maintains approximate orthogonality to support sparse density-kernel calculations. It uses FOE for the very large systems considered.

B. Optimization of the support functions and density kernel

  • The target function combines band-structure energy minimization for accuracy with confinement for support-function localization.The confinement prefactor is smoothly reduced during the run, leaving band-structure energy minimized at the end.Source: B. Optimization of the support… · p. 3
  • Support functions are kept approximately orthogonal because this makes the overlap matrix close to identity and helps preserve sparsity in the density kernel.Orthogonality makes dual and nondual support functions identical, allowing their cutoffs to be chosen along similar lines.Source: B. Optimization of the support… · p. 3
  • Quasi-orthogonal support functions help retain sparsity in the matrix product used by FOE to determine the density kernel.The inverse square root of the overlap matrix is cheaper to calculate by polynomial expansions when the overlap matrix is close to identity.Source: B. Optimization of the support… · pp. 3–4
  • Approximate rather than strict orthogonality avoids potential convergence problems because orthogonality and locality generally conflict.A Lagrange multiplier approximately preserves orthogonality after initial explicit orthogonalization; off-diagonal overlap elements are typically about 0.1.Source: B. Optimization of the support… · p. 4
  • The paper uses FOE to determine the density kernel because it is the offered method that allows calculations for very large systems.Source: B. Optimization of the support… · p. 4

III. ASSESSMENT OF THE ACCURACY

The accuracy assessment compares the linear-scaling code with the cubic version of BigDFT as a reference. Calculations use free boundary conditions and PBE unless otherwise stated.

III. ASSESSMENT OF THE ACCURACY

  • The cubic version of the code7 is used as the reference for assessing the linear-scaling version's precision.Source: §III · p. 4
  • The paper treats reproduction of the cubic version's results as evidence that the linear version has reached high precision.The cubic version is described as having been shown to yield very accurate results.Source: §III · p. 4
  • Calculations use free boundary conditions and the PBE functional unless otherwise stated.Source: §III · p. 4

A. Comparison of energies and forces

Across six systems, the linear and cubic versions show a roughly constant energy offset and small force-norm differences. Similar input parameters are used across systems, with support-function cutoffs varied by element.

A. Comparison of energies and forces

  • The six-system comparison uses similar parameters, with support-function cutoff radii as the main system-dependent setting.Cutoffs are selected by periodic-table row, and the kernel cutoff is set to 9.0 bohr.Source: A. Comparison of energies and forces · p. 4
  • The chosen parameters are described as fairly universal, allowing other systems to be calculated without adjustments.Source: A. Comparison of energies and forces · p. 4
  • About 10 meV/atom energy offset and typically 10^-3 hartree/bohr force-norm difference separate linear from cubic BigDFT in the six-system comparison.These differences were obtained without fine-tuning; support-function optimization also removes Pulay forces, making force evaluation straightforward.Source: A. Comparison of energies and forces · p. 4, Figure 2 · p. 4

B. Calculations of energy differences

For the 660-atom silicon wires, the linear method reproduces the pure-to-impure energy difference much more closely than its absolute energies. The mismatch is 0.15% of the reference difference, as the absolute-energy offset cancels to a large degree.

B. Calculations of energy differences

C. Geometry optimizations

In a 302-atom alkane geometry optimization, the linear and cubic versions follow essentially identical trajectories. Their energy offset remains nearly constant, and their forces differ by only about 10^-4 hartree/bohr.

C. Geometry optimizations

  • The test used a 302-atom alkane starting from a non-equilibrated structure, with the linear and cubic versions optimized from the same geometry.Source: C. Geometry optimizations · p. 5
  • The optimization's force plateau reflects bond-length adjustments beginning at both alkane endpoints and propagating toward the center.Source: Figure 4 · p. 6
  • 10^-4 hartree/bohr is the force variation between the linear and cubic versions, whose geometry-optimization trajectories are identical.The energy offset remains about 10 meV/atom throughout the run.Source: C. Geometry optimizations · pp. 5–6, Figure 4 · p. 6

D. Energetic ordering

For B12C48, finite-cutoff linear scaling preserves the ground state but reorders excited states, including giving the wrong gap to the first local minimum. This delicate ordering is also sensitive to the functional used.

D. Energetic ordering

  • The approach is not universally applicable when finite-cutoff noise exceeds the energy signal; results then require careful interpretation.Source: D. Energetic ordering · p. 6
  • The linear version retains the B12C48 ground state but reorders excited states and gives the wrong gap to the first local minimum.The ordering requires resolving differences of about 1 meV/atom, which the cubic version can capture.Source: D. Energetic ordering · p. 6, Figure 5 · p. 6
  • Changing the cubic calculation from PBE to LDA leaves the ground and first-excited levels identical but completely jumbles higher levels.Source: D. Energetic ordering · p. 6

IV. SCALING WITH SYSTEM SIZE

The linear-scaling version shows linear runtime and memory growth, while the cubic version rises more steeply; CPU minutes per atom level off beyond a system-size threshold. The authors report moderate resource requirements for large single-point calculations.

  • The linear version exhibits strict linear scaling and is faster even for the smallest water droplet, while the cubic version does not permit calculations beyond about 2000 atoms.The estimated runtime crossover is about 200 atoms; the benchmark covers water droplets from 600 to 15000 atoms.Source: §IV · pp. 6–7, Figure 7 · p. 7
  • A complete single-point calculation for 10000 atoms requires only a few thousand CPU hours in the reported benchmark.The authors characterize this resource use as moderate.Source: §IV · p. 7
  • CPU minutes per atom remain constant above a critical size for DNA fragments and water droplets, with higher values for water droplets because of their geometry.The figure reports CPU minutes per atom for fully self-consistent calculations.Source: Figure 8 · p. 7, §IV · p. 8

V. PARALLELIZATION

For a 14300-atom DNA fragment, parallel speedup is good up to a few thousand cores, with diminishing returns at higher core counts. Memory scaling remains close to ideal, and the best balance depends on system size and computational resources.

A. Parallel efficiency

  • The code uses mixed MPI and OpenMP parallelization; fewer MPI tasks can improve load balancing, reduce communication overhead, and lower memory peaks per compute node.The authors note these benefits can matter when memory is a critical resource.Source: A. Parallel efficiency · p. 8
  • CPU minutes per atom indicate very good speedup up to a few thousand cores, but increase at larger core counts as computation per core becomes too small relative to communication.Load balancing also becomes more difficult in this range, leaving many cores idle much of the time.Source: A. Parallel efficiency · p. 8
  • The high-core-count speedup degradation shifts toward larger core counts as system size increases.The authors attribute the degradation to the small number of operations per core and communication overhead.Source: A. Parallel efficiency · p. 8
  • For the 14300-atom DNA fragment, memory parallelization comes close to perfect scaling, and total walltime continues to decrease even at the largest core count.The runs span 1280 to 25600 cores, using eight OpenMP threads per MPI task.Source: A. Parallel efficiency · p. 8, Figure 9 · p. 8

B. Parallelization of the sparse matrices

Matrix taskgroups replace wasteful global reductions with reductions over the matrix portions needed by groups of MPI tasks. The reported timings show substantial speedups for matrix construction and density-kernel assembly.

B. Parallelization of the sparse matrices

  • The taskgroup approach avoids a global reduction by restricting reductions to MPI taskgroups holding the matrix portions their members need.Each MPI task belongs to at most two taskgroups, reducing the matrix data each task must hold.Source: B. Parallelization of the sparse… · p. 9
  • The approach addresses a runtime and memory bottleneck in the transposed matrix-construction approach, where a global reduction is unnecessary because each MPI task needs only a small matrix portion.Non-blocking collective MPI routines allow a task to participate in two reductions without serializing the code.Source: B. Parallelization of the sparse… · p. 9
  • Matrix taskgroups accelerate overlap and Hamiltonian matrix calculation by more than a factor of 3, and density-kernel assembly by a factor of 16.The table reports speedups for these two operations; the eight-task illustration shows the taskgroup concept.Source: B. Parallelization of the sparse… · p. 9, Table III · p. 9, Figure 10 · p. 9

VI. CONCLUSION

The authors conclude that localized, quasi-orthogonal support functions adapted in situ enable accurate DFT calculations with strict linear scaling. They report reduced computational resource requirements and accuracy without tedious parameter fine-tuning.

  • Strict linear scaling with localized, quasi-orthogonal support functions adapted in situ is achieved without loss of precision from the cubic approach.The support functions use an underlying Daubechies wavelet basis with compact support.Source: §VI · p. 9
  • The approach reduces the degrees of freedom needed for accurate calculations, lowering the scaling prefactor and the required computational resources.The authors say standard input parameters avoid the tedious fine-tuning often needed by other O(N) approaches.Source: §VI · p. 9
  • A complete single-point calculation for a 10000-atom system requires on the order of one thousand CPU hours, and CPU minutes per atom offers a way to estimate resources for systems of any size.The authors also connect low resource consumption and efficient parallelization to more expensive tasks, though the passage ends before specifying them.Source: §VI · p. 9
Loading 1501.05884v1…