We investigate the source of error in the Thomas–Fermi–von Weizsäcker (TFW) density functional relative to Kohn–Sham density functional theory (DFT). In particular, through numerical studies on a range of materials, for a variety of crystal structures subject to strain and atomic displacements, we find that while the ground state electron density in TFW orbital-free DFT is close to the Kohn–Sham density, the corresponding energy deviates significantly from the Kohn–Sham value. We show that these differences are a consequence of the poor representation of the linear response within the TFW approximation for the electronic kinetic energy, confirming conjectures in the literature. In so doing, we find that the energy computed from a non-self-consistent Kohn–Sham calculation using the TFW electronic ground state density is in very good agreement with that obtained from the fully self-consistent Kohn–Sham solution.

## I. INTRODUCTION

^{1,2}is one of the most widely used

*ab initio*methods in physical, chemical, and material science research for understanding and predicting the properties of material systems. Its conceptual foundation lies in the Hohenberg–Kohn (HK) theorem,

^{3}which states that the total energy of the system is a unique, albeit unknown, functional of its density. This rather formal mathematical concept was made practically useful by the Kohn–Sham (KS) formalism,

^{4}wherein the real system of interacting electrons is replaced by a fictitious system of non-interacting fermions that generate the same electronic density. In particular, the electronic kinetic energy is no longer an explicit functional of the density, but rather takes the form (in atomic units)

^{1}This severely restricts the range of systems that can be studied using KS-DFT.

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.^{5} 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}.

*T*

_{s}was suggested long before the advent of DFT by Thomas

^{6}and Fermi:

^{7}

^{,}

*λ*= 1 was derived by von Weizsäcker.

^{8}This term by itself represents a lower bound on

*T*

_{s}, being exact in the limit of small and rapid (large wavevector) density variations.

^{9}A similar formula with weight factor

*λ*= 1/9 was derived by Kirzhnits,

^{10}the kinetic energy

*T*

_{s}=

*T*

_{TF}+

*T*

_{λW}being exact in the opposite limit of slow, but not necessarily small variations. Realizing the importance of satisfying the constraint arising from the aforementioned lower bound, a number of other local kinetic energy functionals have recently been proposed, including those of the generalized gradient approximation (GGA)

^{11–13}and the Laplacian-level

*meta*-GGA

^{14,15}varieties. Though these functionals have found significant success, they are semi-empirical, in that they involve parameters that could be material dependent. There have also been suggestions to use alternate values for the weight factor,

^{2}which can be made position dependent;

^{16}however, such strategies also depend upon the class of system under consideration and so lack universality.

*T*

_{λW}signals the inability of the local Thomas–Fermi–von Weizsäcker (TFW) functional, i.e.,

*T*

_{s}=

*T*

_{TF}+

*T*

_{λW}, to describe the kinetic energy even for small density variations, if the scale of variations is neither particularly large nor particularly small. This led to thinking, as early as three decades ago,

^{17}of the need for functionals that are nonlocal in the coordinate space, i.e., they depend on the density correlation at finite distances, e.g.,

*K*(

*ρ*(

**r**),

*ρ*(

**r**′)) is selected in such a way as to improve the description of the noninteracting susceptibility $\chi 0\u22121(r,r\u2032)=\delta 2Ts/\delta \rho (r)\delta \rho (r\u2032)$ to reproduce the uniform electron gas limit given by the analytic Lindhard formula.

^{18}Note that a slightly different form, essentially equivalent to that above, has also been proposed:

^{17}

^{,}$Ts=TTF+TW[\rho \u0303(r)]$, where $\rho \u0303(r)=\u222b\rho (r\u2032)K(\rho (r),\rho (r\u2032))dr\u2032$. In subsequent work,

^{19}it was pointed out that a form more consistent with the idea of generalizing the von Weizsäcker–Kirzhnits term to the domain of nonlocal functionals should include log-density gradients in powers adding up to 2, e.g.,

^{20}who represented the kinetic energy as

*a*+

*b*= 5/3. This was later generalized to density-dependent kernels by Wang, Govind, and Carter (WGC):

^{21}

^{,}

*a*+

*b*= 8/3. While showing good results for particular problems,

^{21–28}such nonlocal functionals have found rather limited use in practice due to greater computational expense and the need for specialized kernels to be developed for different material systems.

^{23,25–27,29}Further advances require a more fundamental understanding of the linear response, and in particular, whether it is the deciding factor in determining the error associated with OF-DFT, as conjectured in previous studies.

^{17,19,20,30}

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,^{31} the well-studied^{32,33} Harris–Foulkes functional,^{34,35} as in the shell correction method,^{36} 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^{39} 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}

## II. SYSTEMS AND METHODS

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^{40} 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 formwhere the perturbation parameter −0.3 ≤(1)$G=1+g0001+g0001+g,$
*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 formwhere the perturbation parameter −0.1 ≤(2)$G=(1+3g)\u2212131+gggg1+gggg1+g,$
*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 formwhere the perturbation parameter −0.1 ≤(3)$G=(1+g)\u221212000(1+g)\u2212120001+g,$
*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.^{43} 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).^{40} In the OF-DFT calculations, we choose the TFW kinetic energy functional with weight factor *λ* = 1/8.^{49} 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^{−4} 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^{−3} 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^{−6} ha/atom.

## III. RESULTS AND DISCUSSION

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.

*E*

_{OF→KS},

*E*

_{KS→OF}, and

*E*

_{OF→OF}, the error being defined with respect to

*E*

_{KS→KS}, we define the following root-mean-square measure, referred to as the Δ-value:

^{50}

*E*∈ {

*E*

_{OF→KS},

*E*

_{KS→OF},

*E*

_{OF→OF}}, and the corresponding Δ ∈ {Δ

_{OF→KS}, Δ

_{KS→OF}, Δ

_{OF→OF}}. Note that the difference in energies from

*g*= 0 is used in the definition of the error, since the reference energy within KS-DFT and OF-DFT is different, a consequence of the different energy functionals, i.e., only differences in energy within the same level of theory are meaningful. We emphasize that we intentionally gauge the OF functional against the KS functional, and not against the experiment.

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).

. | Δ_{OF→KS}
. | Δ_{KS→OF}
. | Δ_{OF→OF}
. |
---|---|---|---|

Volumetric strain | 0.001 024 3 | 0.009 120 2 | 0.016 636 8 |

Rhombohedral strain | 0.000 046 6 | 0.002 842 0 | 0.004 405 3 |

[001] Uniaxial strain | 0.000 042 3 | 0.002 093 5 | 0.003 697 6 |

Atomic perturbation | 0.000 024 6 | 0.000 159 8 | 0.000 220 4 |

. | Δ_{OF→KS}
. | Δ_{KS→OF}
. | Δ_{OF→OF}
. |
---|---|---|---|

Volumetric strain | 0.001 024 3 | 0.009 120 2 | 0.016 636 8 |

Rhombohedral strain | 0.000 046 6 | 0.002 842 0 | 0.004 405 3 |

[001] Uniaxial strain | 0.000 042 3 | 0.002 093 5 | 0.003 697 6 |

Atomic perturbation | 0.000 024 6 | 0.000 159 8 | 0.000 220 4 |

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.

^{1,51}

^{3}Indeed, in the case of strains and atomic perturbations, the changes in energy are dictated by the elastic constants and interatomic force constant matrix, respectively, which can be derived in the context of the above formalism, as done previously in KS-DFT.

^{51}

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,^{52} 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.^{28}

## IV. CONCLUDING REMARKS

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

## REFERENCES

*Electronic Structure: Basic Theory and Practical Methods*

*Horizons of Quantum Chemistry*

*Handbook of Materials Modeling*

*Mathematical Proceedings of the Cambridge Philosophical Society*

*Principles of the Theory of Solids*

We have also tested other choices of *λ* and have found that the final conclusions remain unchanged