Assessing the source of error in the Thomas–Fermi–von Weizsäcker density functional

An alternative to replacing the system of interacting electrons with a fictitious system of non-interacting fermions is to replace it instead with a fictitious system of non-interacting bosons. This can be achieved by approximating the kinetic energy T_s using an explicit functional of the density, the resulting formalism referred to as orbital-free (OF) DFT.⁵ This generally amounts to solving a Schrödinger-type equation for only one electronic state, which corresponds to the square root of the density. Though the computational cost of OF-DFT scales linearly with system size, thus overcoming the cubic-scaling bottleneck of KS-DFT, it is rarely used in practice due to the lack of accurate approximations for T_s.

An alternative strategy to nonlocal density functionals is to perform a non-self-consistent KS calculation with the electronic ground state density from OF-DFT as input. In particular, the ground state electron density obtained from an OF-DFT calculation is used as input to a single self-consistent field (SCF) iteration of a KS calculation, the quantities so obtained then being used to compute the energy, either using the Kohn–Sham functional,³¹ the well-studied^32,33 Harris–Foulkes functional,^34,35 as in the shell correction method,³⁶ or a weighted combination of the two, as in the orbital-corrected orbital-free (OO) DFT.^37,38 The formulations/implementations^31,36,38 related to the Strutinsky shell correction method³⁹ are for molecular systems in the context of TF OF-DFT, while the focus in the current work is on TFW and condensed matter systems. The OO-DFT method, developed in the context of WGC OF-DFT, displays very good agreement with KS-DFT. Though interesting, it was, however, semi-empirical, with dependence on a parameter that determines the relative contributions of the Kohn–Sham and Harris–Foulkes energies, in contrast to the above-mentioned nonlocal functionals underpinned by the linear response theory. Furthermore, WGC already incorporates some amount of linear response, complicating the ability to draw clear inferences. Finally, the numerical evidence consists of just two cases, namely, the energy–volume curves for fcc Ag and cubic-diamond Si. To what extent this idea is generally applicable, and how (if at all) it is related to the linear response theory, remain to be clarified.

In this paper, we address the aforementioned issues. First, we verify that the “single-shot” KS calculation with TFW OF-DFT ground state density as input is close to the fully self-consistent KS-DFT solution for a range of materials, including different crystal structures subject to volumetric and symmetry-lowering perturbations. Second, we argue that the success of such a strategy indicates that the main limitation of TFW OF-DFT is indeed its inability to properly describe the correct linear response, as conjectured in the literature.^17,19,20,30

We consider body-centered cubic (BCC), face-centered cubic (FCC), hexagonal close packed (HCP), and body-centered tetragonal (BCT) crystals of magnesium (Mg), aluminum (Al), and indium (In), as well as diamond cubic (DC) and hexagonal diamond (DH) (also known as lonsdaleite) crystals of silicon (Si). These systems form a diverse set that includes a simple metal, a transition metal, and a semiconductor in a variety of lattice configurations. Importantly, well-tested local pseudopotentials⁴⁰ are available for the chemical elements in question, i.e., Mg, Al, In, and Si, allowing for a careful comparison of the results obtained from KS-DFT and OF-DFT calculations.

Unless specified otherwise, we choose the primitive unit cells for each of the systems, i.e., one-atom cells for the BCC and FCC lattices, two-atom cells for the HCP, BCT, and DC lattices, and four-atom cells for the DH lattice. For the HCP and BCT lattices, we choose the ideal c/a ratios 1.633 and 1.414, respectively. After determining the equilibrium configurations, we consider the following strains and atomic displacements:

• Volumetric strains for each of the aforementioned systems. In this case, the strain tensor takes the form

  $G=1+g0001+g0001+g,$

  (1)

  where the perturbation parameter −0.3 ≤ g ≤ 0.3, with g < 0 and g > 0 corresponding to the contraction and expansion of the unit cell, respectively.
• Symmetry-lowering, volume-preserving, rhombohedral strains for Mg, Al, and In in the FCC and BCC crystal configurations and Si in the DC crystal configuration. In this case, the strain tensor takes the form

  $G=(1+3g)−131+gggg1+gggg1+g,$

  (2)

  where the perturbation parameter −0.1 ≤ g ≤ 0.1, with g < 0 and g > 0 corresponding to the compression and elongation of the unit cell, respectively, along the $111$ direction.
• Volume-preserving uniaxial strains along the $001$ direction for Mg, Al, and In in the HCP and BCT crystal configurations and Si in the DH crystal configuration. In this case, the strain tensor takes the form

  $G=(1+g)−12000(1+g)−120001+g,$

  (3)

  where the perturbation parameter −0.1 ≤ g ≤ 0.1, with g < 0 and g > 0 corresponding to the compression and elongation of the unit cell, respectively.

• Symmetry-lowering atomic perturbations (i.e., frozen phonons) for Al and Mg in the BCC, FCC, and HCP crystal configurations and In in the FCC, HCP, and BCT crystal configurations. For BCC and FCC, we choose the two-atom conventional cell and two-atom tetragonal cells, respectively, rather than the one-atom primitive cell used in other simulations. In all cases, the atom is perturbed along the $001$ direction, with the z coordinate of one of the atoms changed as z → (1 + g)z, where the perturbation parameter −0.05 ≤ g ≤ 0.05.

The systems considered here ensure that both slow and rapid variations in the density are encountered, especially when the internal parameters are varied, enabling a thorough analysis of the error in the TFW functional. Note that even in the case of uniform strains, there are rapid changes in the electron density for many of the systems, as evidenced by the large errors in the TFW OF-DFT energies (supplementary material).

All calculations are performed using the M-SPARC code,^41,42 which is a MATLAB version of the large-scale parallel electronic structure code, SPARC.⁴³ It employs the real-space finite-difference method, whose formulation and implementation in the context of KS-DFT and OF-DFT can be found in previous studies.^44–47 We employ the local density approximation (LDA)^4,48 for the exchange-correlation functional and use the bulk-derived local pseudopotentials (BLPs).⁴⁰ In the OF-DFT calculations, we choose the TFW kinetic energy functional with weight factor λ = 1/8.⁴⁹ In the KS-DFT calculations, we perform Brillouin zone integration using a 15 × 15 × 15 Monkhorst–Pack grid for the FCC, BCC, DC, and DH lattices, and a 15 × 15 × 10 grid for the HCP and BCT lattices, which ensures that the energies are converged to within 10⁻⁴ ha/atom. In all calculations, we employ a 12th order finite-difference approximation and a grid spacing of 0.4 bohr, which ensures that the computed energies are converged to within 10⁻³ ha/atom. Finally, the change in energy arising due to a perturbation, which is the main quantity of interest in the present work (Sec. III), is converged to within 10⁻⁶ ha/atom.

We use the framework described in Sec. II to perform KS-DFT and OF-DFT calculations for the selected systems. In particular, for each system, we compute the four energies listed below:

• E_KS→KS ≔ E_KS(ρ_KS): KS-DFT energy E_KS corresponding to the KS-DFT ground state density ρ_KS. This is obtained by performing a standard electronic ground state calculation in KS-DFT.

• E_OF→KS ≔ E_KS(ρ_OF): KS-DFT energy E_KS corresponding to the OF-DFT ground state density ρ_OF. This involves the calculation of orbitals for the given ρ_OF, i.e., a single self-consistent field (SCF) iteration in KS-DFT.

• E_KS→OF ≔ E_OF(ρ_KS): OF-DFT energy E_OF corresponding to the KS-DFT ground state density ρ_KS.

• E_OF→OF ≔ E_OF(ρ_OF): OF-DFT energy E_OF corresponding to the OF-DFT ground state density ρ_OF. This is obtained by performing a standard electronic ground state calculation in OF-DFT.

In Fig. 1 and Table I, we summarize the results so obtained, with the detailed data available in the supplementary material. We observe the following trend in the Δ-errors: Δ_OF→KS < Δ_KS→OF < Δ_OF→OF. In particular, the values of Δ_OF→KS are quite small, which suggests that the ground state density in OF-DFT is close to that in KS-DFT. Furthermore, since the values of Δ_KS→OF are relatively large, it can be inferred that the energy errors in OF-DFT are not a consequence of the errors in the ground state density, but rather due to a fundamental limitation in the energy functional. Indeed, similar trends and inferences follow when comparing physical observables such as equilibrium volume and bulk modulus (supplementary material).

Though the current findings cannot be considered universal, we expect them to be generally true for TFW OF-DFT, given that the chosen systems contain several families of materials (metal and semiconductors), classes of structures (close-packed, open, and intermediate), symmetries (cubic, hexagonal, tetragonal, orthorhombic, and rhombohedral), and types of distortions (cell and atom). To further bolster our findings, we look to theory, as described below.

Since the difference between ΔE_KS→KS and ΔE_OF→KS is small, as found in the numerical results above, it follows from Eqs. (5) and (6) that the ground state densities of KS-DFT and OF-DFT are close, which also justifies the use of the linear response theory. Furthermore, since the difference between ΔE_KS→OF and ΔE_KS→KS is relatively large, as found in the numerical results above, it follows from Eqs. (5) and (7) that the error in the energy for OF-DFT calculations is due to the poor representation of the linear response susceptibility χ_OF relative to χ_KS, which confirms the need for developing alternate functionals with better linear response. Indeed, as discussed in Sec. I, this has motivated the development of a number of local^11–15 and nonlocal^{17,19–21,23,25–27,29} functionals. Given that these functionals are generally more accurate than TFW,⁵² it is expected that the current findings are also applicable to them, though the accuracy of the ground state density and computed energy may vary between different functionals. Since the focus of the current work is to investigate the source of error in the TFW functional, rather than provide a comparison between different functionals, we refrain from such comparison here.

The above results also suggest a possible strategy to accelerate KS-DFT calculations without significant loss of accuracy. In particular, the ground state density computed from OF-DFT, which scales linearly with system size, can be used as input to accelerate SCF convergence or just perform a single SCF iteration in KS-DFT. Indeed, to increase the fidelity of such single-shot calculations, nonlocal pseudopotentials, which are the standard in KS-DFT, need to be incorporated into OF-DFT.²⁸ 

In this work, we have systematically investigated the source of error arising in the TFW density functional relative to the Kohn–Sham DFT. In particular, through numerical studies on a variety of materials, for a range of crystal structures subject to strains and atomic displacements, we have found that while the ground state electron density in the TFW variant of orbital-free DFT is close to the Kohn–Sham ground state density, the corresponding energy differs significantly from the Kohn–Sham value. We have shown that these differences arise due to the poor representation of the linear response within the TFW approximation for the electronic kinetic energy, therefore confirming conjectures in the literature. In so doing, we have found that the energy computed from a non-self-consistent Kohn–Sham calculation using the TFW ground state density as input is in very good agreement with the energy obtained from the fully self-consistent Kohn–Sham solution.

The development of more general and accurate electronic kinetic energy functionals for use in orbital-free DFT, possibly aided by state-of-the-art machine learning techniques, is therefore a worthy subject of future work.

Plots of variation in ΔEKS→KS, ΔEOF→KS, ΔEKS→OF, and ΔEOF→OF with perturbation.

The authors would like to thank Sam Trickey for helpful discussions on orbital-free DFT literature. The authors also sincerely acknowledge the invaluable help from Maria Emelyanenko. The authors also thank the anonymous reviewers for their valuable comments and suggestions. B.T. acknowledges the support of the Quantum Science and Engineering Center (QSEC) at George Mason University. X.J., P.S., and J.P. acknowledge the support of Grant No. DE-SC0019410, funded by the U.S. Department of Energy, Office of Science. This work was performed, in part, under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory, under Contract No. DE-AC52-07NA27344.

The authors have no conflicts to disclose.

Bishal Thapa: Data curation (lead); Investigation (equal); Methodology (supporting); Software (supporting); Validation (equal); Writing – original draft (lead). Xin Jing: Methodology (equal); Software (lead). John E. Pask: Funding acquisition (supporting); Supervision (supporting); Writing – review & editing (supporting). Phanish Suryanarayana: Data curation (supporting); Formal analysis (equal); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Resources (supporting); Software (supporting); Supervision (supporting); Writing – original draft (supporting); Writing – review & editing (lead). Igor I. Mazin: Conceptualization (lead); Data curation (supporting); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Resources (lead); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (supporting).

The data that support the findings of this study are available within the article and its supplementary material.