Imposing Correct Jellium Response Is Key to Predict the Density Response by Orbital-Free DFT

Orbital-free density functional theory (OF-DFT) constitutes a computationally highly effective tool for modeling electronic structures of systems ranging from room-temperature materials to warm dense matter. Its accuracy critically depends on the employed kinetic energy (KE) density functional, which has to be supplied as an external input. In this work we consider several nonlocal and Laplacian-level KE functionals and use an external harmonic perturbation to compute the static density response at T=0 K in the linear and beyond linear response regimes. We test for the satisfaction of exact conditions in the limit of uniform densities and for how approximate KE functionals reproduce the density response of realistic materials (e.g., Al and Si) against the Kohn-Sham DFT reference which employs the exact KE. The results illustrate that several functionals violate exact conditions in the UEG limit. We find a strong correlation between the accuracy of the KE functionals in the UEG limit and in the strongly inhomogeneous case. This empirically demonstrates the importance of imposing the limit of UEG response for uniform densities and validates the use of the Lindhard function in the formulation of kernels for nonlocal functionals. This conclusion is substantiated by additional calculations for bulk Aluminum (Al) with a face-centered cubic (fcc) lattice and Silicon (Si) with an fcc lattice, body-centered cubic (bcc) lattice, and semiconducting crystal diamond (cd) state. The analysis of fcc Al, and fcc as well as bcc Si data follows closely the conclusions drawn for the UEG, allowing us to extend our conclusions to realistic systems that are subject to density inhomogeneities induced by ions.


I. INTRODUCTION
First-principles methods based on the electronic density are often used for simulations of electronic structures in physics and chemistry.Usually, there is a punishing correlation between the increase in the accuracy of a method and its computational cost.Indeed, there is a computational bottleneck for methods such as Kohn-Sham density functional theory (KS-DFT) or, more accurately, generalized KS-DFT with hybrid exchangecorrelation (XC) functionals, which strongly hinders the simulation of large systems (e.g., with the number of particles of the order of N ∼ 10 4 and more).The simulation of large systems is considered to be important for the adequate description of processes involving large numbers of particles such as phase transitions [1], nucleation [2], freezing or melting [3,4], collective ion oscillation (normal) modes [5][6][7], and to calculate transport properties such as diffusion and viscosity [8][9][10].These properties are often accessed by combining molecular dynamics (MD) simulations of ions with the electron forces computed using a density functional theory based method.
Orbital-free density functional theory (OF-DFT) [11] is one of the DFT approaches being actively developed * z.moldabekov@hzdr.de† xuecheng.shao@rutgers.edu‡ m.pavanello@rutgers.edufor the simulation of large systems because of its generally linearly scaling computational cost with respect to the number of particles.OF-DFT has been employed successfully for the description of materials [11], melted metals [12], and even nanoparticles [13].Additionally, the computational cost of OF-DFT is not sensitive to the variation in temperature since it does not use orbitals.This is an important computational advantage of OF-DFT when it comes to the simulation of phenomena at high temperatures [14], at which thermal KS-DFT simulations require large number of orbitals to correctly capture thermal excitations.Furthermore, being computationally inexpensive, OF-DFT can be used as an intermediate step to accelerate KS-DFT based MD simulations by optimizing an initial ionic configuration [15].Moreover, the OF-DFT approach is one of the main tools for the simulation of near-mesoscopic-scale dynamics taking into account the quantum nonlocality in plasmonics [16,17], quantum plasmas [18][19][20][21] as well as plasmonics in nanomaterials [22][23][24].Recent advances have extended the applicability of OF-DFT to dynamic structure factors [25] and optical properties of clusters and nanostructures [26][27][28].
In contrast to KS-DFT, the non-interacting kinetic energy (KE) functional in OF-DFT is not expressed as a pure functional of the KS orbitals but it is instead expressed as a pure functional of the electron density.The functional is known to exist, but its analytical pure dependency on the electron density is not known exactly arXiv:2304.11169v2[physics.comp-ph]18 Oct 2023 and has to be approximated.Therefore, the key problem is to design a KE functional, and the corresponding non-interacting free energy at finite temperatures.Following early seminal works by, e.g., Thomas and Fermi (TF) [29,30], Kohn and Sham [31], Perrot [32], and Wang and Teter (WT) [33], it became a standard to design KE functionals that reproduce the correct density response function in the limit of the ideal uniform electron gas (UEG) [34][35][36][37][38].This connection to the archetypical freeelectron system delivers a certain level of universality, which, e.g., can lead to accurate results for bulk properties of metals and semiconductors [13,[39][40][41][42][43].
Using the UEG model, a commonly used relation (constraint) for the construction of a KE functional T s [n] is based on the Lindhard function χ Lin : [19,31,44] −F δ 2 T s δn(r)δn(r where η = k/ 2k 0 F is a wavenumber in the units of the Fermi wavenumber k 0 F = (3π 2 n 0 ) 1/3 (with n 0 being the mean electron density) and F[...] denotes the Fourier transformation operator.For example, so-called semilocal KE functionals can be build by using the long wavelength expansion of Eq. (1) at η < 1 [39].
The engineering of advanced KE functionals is accompanied by involved theoretical manipulations with the constrain that the Lindhard function is reproduced if one first takes the UEG limit n → n 0 and then performs the Fourier transform according to the relation (1).For example, using a series expansion of χ −1 Lin (η) for η < 1, one recovers a constant term (associated to the Thomas-Fermi functional) and corrections which come with powers of η = kn0 2k F n0 and η 2 = k 2 n0 4k 2 F n0 whose real-space counterparts are s(r) = |∇n(r)|/ (2k F (r)n(r)) and q(r) = |∇ 2 n(r)|/ 4k 2 F (r)n(r) (where k F (r) = (3π 2 n(r) 1/3 is the local Fermi wavenumber).Another example of the operations connecting a real space density inhomogeneity with the density perturbation wavenumber in χ −1 Lin (η) is the functional integration using the local density to map system properties locally on the UEG at the corresponding density [42,43].These schemes for the engineering of KE functionals based on Eq. ( 1) can be symbolically expressed as Being based on assumptions that can be theoretically justified to a certain degree, the KE functionals constructed on top of the Lindhard function might not be traceable back numerically to Eq. (1) in the limit of the UEG leading to an inconsistency: Additionally, nonlocal KE functionals derived from Eq. (1) will suffer from bias towards describing uniform densities best.It is yet unclear whether these functionals are capable of producing accurate density response for the UEG and realistic materials in the linear and beyondlinear response regimes.
In this work, we test some of the open questions mentioned above by implementing a direct perturbation approach.The electronic system is harmonically perturbed allowing us to compute the static density response function of UEG and any bulk material using OF-DFT as well as KS-DFT.This direct perturbation approach constitutes an effective tool to check inconsistencies in the sense of Eq. ( 3) and generally the quality of the density response in the linear and beyond-linear regimes.We specifically focus on the case of fully degenerate electrons with a temperature T ≪ T F , where T F is the Fermi temperature.
A particular focus of this work is on KE functionals designed imposing Eq. ( 1).Since the Lindhard function based construction is not optimal for bound states and systems with band gaps, other strategies such as utilizing the jellium-with-gap model [45], constrains on the Pauli potential and constraints on atom densities (e.g.Kato cusp condition) have been used.These observations motivate us to employ the method of the direct perturbation approach to compute the static density (linear and beyond-linear) response function for various KE functionals to inspect the severity of the bias introduced by Eq. ( 1).We reiterate that, in addition to the non-interacting kinetic energy (KE) functional, the XC functional needs to be supplied as an input of any DFT simulation.However, because we focus here on KE functionals, we purposely do not test the variability of our results with respect to different choices of XC functionals.
A particularly interesting problem to be analyzed using the direct perturbation method is related to the applicability of the UEG based KE functional when the density perturbation is beyond the linear response regime.Indeed, when the perturbation is strong enough, linear density response theory (like χ Lin ) can be grossly inaccurate for the description of the density perturbation [46][47][48][49].Here we show that one can analyze the performance of the KE functionals by observing of the response from approximate OF-DFT KE functionals against KS-DFT results with increasing amplitude of the harmonic perturbation.
We consider an array of KE functionals, for example the nonlocal WT functional [33] and MGP functionals [40] (the latter uses functional integration techniques to define the kernel of the nonlocal term), the Laplacianlevel, meta-GGA PGSL functional [41], and LWT and LMGP functionals created by introducing a local density dependence into the kernels of MGP and WT, respectively [43].The constraint (1) was used for the construction of these functionals.By design, the WT and LWT (and to an extent also MGP and LMGP) should satisfy Eq. ( 1), in principle, for all wavenumbers and PGSL should satisfy Eq. ( 1) at η ≲ 0.5 (i.e., k ≲ k 0 F ).The application of the direct perturbation approach to real materials is provided for several phases of bulk Alu-minum (Al) and Silicon (Si).We analyzed the density response of the valence electrons to a weak external harmonic field at different wavenumbers.This allowed us to relate the performance of the KE functionals in the UEG limit to the quality of the density response description for these realistic, inhomogeneous materials.
Finally, taking into account the interest in the application of OF-DFT for warm dense matter research [19-21, 25, 50-52], we formulate recommendations for the application of these functionals at such extreme densities and temperatures for degenerate electrons.
The paper is organized as follows: In Sec.II, we provide a brief theoretical background of the considered KE functionals to emphasise the relevant physical aspects.We also introduce the direct perturbation method.We provide the simulation details in Sec.III.The results and discussions are presented in Sec.IV where we present static (linear and beyond-linear) response simulations for the UEG and bulk Al and Si.The paper is concluded by summarizing main findings and providing an outlook over potential future works.

II. THEORETICAL BACKGROUND
In this section, we first briefly introduce the nonlocal and Laplacian-level KE functionals considered in this work.Then, we discuss how the constraint (1) is used in their construction.Lastly, we introduce the harmonic perturbation approach.
Other nonlocal KE functionals based on the Lindhard function and considered in this work are MGP [40], LMGP and LWT [43].In contrast to the WT model for which the KE potential δT WT /δn is computed using an ansatz (4) for the KE functional, the MGP potential is calculated by functional integration [40].The difference between MGP and LMGP is the way how the Lindhard function is used.In the MGP, χ −1 Lin (k) is computed using the mean density n 0 of the valence electrons.In the case of the LMGP, χ −1 Lin (k) is computed for each grid point using the density value n 0 → n(r) on this grid point.Similarly, the LWT potential is computed using an n 0 → n(r) mapping locally for the WT KE kernel (4) [43].Therefore, being nonlocal functionals based on the Lindhard function and with rather non-trivial mathematical manipulations, the MGP, LMGP, and LWT functionals are perfect candidates to showcase the utility of the direct perturbation approach for computing the density response and for testing functionals (i.e. in the limit of the UEG).
The PGSL KE functional (10), by its design, should reproduce the Lindhard density response function of the UEG at k < 2k 0 F .In Ref. [41], it was demonstrated that the PGSL KE functional provides accurate results for the bulk properties of metals and semiconductors without using system-dependent parameters.
The second semilocal KE functional that we use is PGS (standing for "Pauli-Gaussian second order") [41], which has the same functional form as Eq. ( 9) but with the following kernel: The PGS functional is designed to satisfy the secondorder gradient expansion following from the long wave length expansion of the inverse Lindhard function [19,32,37,54].
The third semilocal KE functional that we consider is the generalized gradient approximation by Luo, Karasiev, and Trickey (LKT) [55]: where with a = 1.3.
The LKT functional is designed to satisfy the positivity of the Pauli functional and potential, and was tested successfully on simple metals and semiconductors [55].

B. Harmonic Perturbations
The OF-DFT calculations of the electron density distribution are performed by minimisation of the total energy under the constraint of a constant particle number.The corresponding Lagrangian is given by where W H [n] is the classical Coulomb repulsion between electrons in a mean-field (Hartree) approximation, W ei [n] = v PP (r)n(r)dr is the potential energy due to the electron-ion interaction (v PP is an ionic pseudopotential), V xc [n] is the XC energy, and µ c is a constant defining the chemical potential at a fixed number of electrons N in the simulation cell.Additionally, the Lagrangian (11) has a term corresponding to the external harmonic perturbation with an amplitude A and wavenumber k.
The latter leads to a deviation of the electron density from its mean value, For the UEG, in the linear response regime (i.e., for weak perturbations), the change in the density δn A,k (r) can be computed using the static electronic density response function χ(k), Therefore, the calculation of the electron density perturbation δn A,k (r) due to an external harmonic field and the inversion of Eq. ( 13) allow one to compute the static density response function χ(k).We note that this method can be generalized to inhomogeneous systems and to the dynamic density response function [56][57][58].
Next, if we neglect the XC energy and set V xc = 0 in Eq. ( 11), the OF-DFT calculation of the harmonically perturbed system delivers the screened non-interacting density response function, with the screening effect being due to the presence of the Hartree term in the Lagrangian (11).For this case, the exact analytical solution in thermodynamic limit reads [36]: where χ s is the static Kohn-Sham response which reduces to χ Lind for the UEG and the inverse reduces to a simple "one-over" operation.The resulting response in Eq. ( 14) is referred to as random phase approximation (RPA).
In this way, using the OF-DFT results for the density response function of the UEG in the RPA and Eq. ( 14), we can directly test whether a given KE functional is able to reproduce the Lindhard function at the relevant wavenumbers.
Increasing the amplitude of the external perturbation A eventually leads to a density perturbation outside of the linear response domain.We recall that the KS-DFT method delivers an exact result for a given choice of XC functional (including omitting it completely to obtain the so-called RPA approximation) .Therefore, by performing a comparative analysis of the OF-DFT and KS-DFT results for the density perturbation beyond the linear response regime, we can rigorously assess to which degree the quality of the KE functionals built on the basis of the linear response function of the UEG deteriorates.In general, OF-DFT is used to describe inhomogeneous systems.Therefore, we will also provide examples for its application to a system that contains ions (bulk Al and Si).For context, in prior works, the direct perturbation approach has been used to compute the static density response function and XC kernel of warm dense matter using quantum Monte-Carlo [46,49] and thermal KS-DFT methods [48,56,[59][60][61].Here, we extend its use to probe the quality of KE density functionals in OF-DFT.
We note that the response functions considered so far are frequency independent, that is, they describe the density response due to an adiabatic perturbation whereby the electrons have an infinite time to adjust to the applied perturbation.The ability of the KE functional approximations to produce quality adiabatic responses is key to quality time-dependent OF-DFT simulations [22,24,28].In Ref. 28 Jiang et al. argue that, while nonlocal functionals (specifically LMGP) provide good approximations to the adiabatic response, semilocal, GGA functionals also deliver an accurate adiabatic response.The test systems in that work were comprised of clusters of various sizes.By considering systems with periodicity, this work aims at providing valuable information as to whether (semi)local KE functionals live up to the expectations set for in Ref. 28.

III. CALCULATION PARAMETERS
The OF-DFT simulations were performed using the DFTpy code that is based on a plane-wave expansion of the electron density [62].The KS-DFT calculations were performed using the ABINIT package [63][64][65][66][67][68].We simulated a free electron gas in the ground state at a characteristic metallic density n 0 = 2 × 10 23 cm −3 with periodic boundary conditions.
For the calculation of the linear response function, the amplitude of the perturbation is set to A = 0.01 (in Hartree), which creates a weak density perturbation described accurately by linear response theory.To test the KE functionals beyond the linear response regime, we consider A = 0.1, A = 0.5, and A = 1.The wavenumber of the external perturbation is set along the x-axis.The wavelength of the perturbation has to be commensurate with the size L of the main simulation cell, which is defined by the relation N L 3 = n 0 .Accordingly, the wavenumber values of the external harmonic perturbation are given by k = j × k min , where k min = 2π/L and j is a positive integer number.We reconstruct the density response function dependence on k in the range 0.1 ≤ k/k 0 F < 4 by performing calculations for different j values and for different numbers of particles.
The OF-DFT calculations of free electron gas were performed for 7168, 66, 38, and 14 electrons in the simulation cell.For 7168 particles we have k min /k 0 F ≃ 0.105.The grid spacing was set to L/200.
The KS-DFT simulations were performed for N = 38 electrons in the main cell with 40 bands and with 30 Ha energy cutoff.The k-point sampling was set to 10 × 10 × 10.We present the results from KS-DFT simulations with zero XC functional and with LDA [69].We crosschecked the KS-DFT results by reproducing our ABINIT data with independent calculations based on the GPAW code [70][71][72][73].
For the simulation of the valence electrons in a bulk of Al and Si, we used the unit primitive cells of the corresponding lattice structures.The presented KS-DFT simulation results are obtained using Quantum ESPRESSO (QE) [74,75] implemented in the Python engine QEpy [76].The energy cutoff was set to 30 Ha (corresponding to a wavefunction plane wave cut-off energy of 816 eV) and the k-point sampling was set to 16 × 16 × 16.The QE results for a considered setup are cross checked by the GPAW simulations.OF-DFT simulations for Al using DFTpy were run with a kinetic energy cutoff of 240 eV.Both KS-DFT and OF-DFT calculations for Al were run with BLPS local pseudopotentials [77].

IV. RESULTS
We first present results for the static linear response function of the electron gas computed using an external harmonic field with amplitude A = 0.01 at wavenumbers 0.1 ≤ k/k 0 F < 4.Then, we present an analysis of the accuracy of OF-DFT calculations for the strongly perturbed (beyond linear response) electron gas.

A. Linear density response
Let us first consider the ideal electron gas neglecting all XC effects.In Fig. 1, we present the results for the density response function computed using different KE functionals and setting the XC functional to zero.In Fig. 1a), we compare the OF-DFT results with the screened (RPA) Lindhard density response function (14).Additionally, in Fig. 1a), we show the data computed using KS-DFT with zero XC functional, but non-zero Hartree mean-field potential.We observe that, as expected, the KS-DFT data for ideal electron gas accurately reproduce the analytical RPA result.The same is the case for the OF-DFT data based on the WT KE functional.In contrast, we see that the MGP KE functional based OF-DFT data significantly deviate from the exact RPA data at 1 ≲ k/k 0 F ≲ 2. The MGP based results show substantial disagreements with the exact RPA result at 0.5 ≲ k/k 0 F ≲ 2. The design choices leading to the MGP kernel are such that the kernel should be optimized (i.e., optimal values for two parameters defining the so-called "kinetic electron" [40]).The results presented here show that without optimization MGP fails to deliver correct response across a wide array of wavevectors.Fig. 1 shows that the LMGP results are in good agreement with the exact RPA data at all considered wavenumbers.This is in line with the design choices that led to the formulation of LMGP which were to yield a KE functional based on MGP with no adjustable parameters [43].Among the nonlocal functionals, we note that the LWT functional (by design) is computed by using the local density approximation in the KE kernel n 0 → n(r) [43].The static density response function computed using 38 and 24 particles show close agreement of the LWT based data with WT based re-sultsIn the case of other considered smaller and larger particle numbers, the LWT became numerically unstable at intermediate wavenumbers (0.5 ≲ k/k 0 F ≲ 2) showing large deviations from the WT based data for the UEG.To avoid possible spurious results, we will not consider LWT for the density response in more non-trivial cases of strong perturbations and real materials.
The PGSL and PGS functional based results are in close agreement with the exact RPA data at k ≲ 2k 0 F and at k ≲ k 0 F .This is in agreement with the analytically enforced constraints to the PGSL and PGS.In general, all considered KE functionals in Fig. 1a) closely agree with the exact solution (14) at k < k 0 F .This can be understood by observing that the screen- ing leads to χ RPA (k) ≈ −k 2 /4π in Eq. ( 14) at small wavenumbers since χ Lin (k → 0) → const [78].Therefore, the Coulomb term 4π/k 2 in the denominator of Eq. ( 14)-representing screening effect-dominates the k dependence of the χ RPA (k) at small wavenumbers.This can mask the inaccuracies of the OF-DFT results at k < k 0 F .In order to assess the quality of the ideal density response function without screening, we invert Eq. ( 14) with χ RPA (k) = χ DFT RPA (k) being computed using the considered KE functionals in OF-DFT and using KS-DFT data.The results obtained in this way are presented in Fig. 1b).Additionally, in Fig. 1b), we compare the DFT results with the Lindhard function which defines the constraint (1) used to construct the considered KE functionals.In the k → 0 limit, all considered KE functionals reproduce the corresponding limit of the Lindhard function as it is illustrated in Fig. 1b) for the point at k = 0.1k 0 F computed using N = 7168 particles.In fact, it is the limit of the TF model.At large wavenumbers k > 2.5k 0 F dominated by the vW KE term, the LWT, MGP, and LMGP KE functionals are able to correctly describe the density response function of the ideal electron gas.This is indeed the well-known single particle regime in which the vW functional provides the exact kinetic energy value.The WT KE functional accurately reproduces the Lindhard function at all considered k values.Among other constrains, the MGP functional was designed to satisfy the relation (1) for, in principle, all wavenumbers.However, from Fig. 1, we see that the MGP results, for the reasons mentioned before, do not satisfy the relation (1) at 0.1 ≲ k/k 0 F ≲ 2. In contrast and in line with the discussion above, the LMGP based ideal density response function is in good agreement with the Lindhard function at the considered wavenumbers.From Fig. 1, and following the expected behavior, the PGSL functional reproduces the Lindhard function for k ≲ 2k 0 F and approaches the Lindhard function at large wavenumbers k > 3k 0 F .In contrast, the PGS functional does not give the correct large wavenumber limit defined by the vW KE functional.Finally, we find close agreement between the LWT and the WT based results.
Next, we consider an interacting electron gas using the LDA XC functional and the corresponding density response function χ LDA (k).In Fig. 2, we show the results for χ LDA (k) computed using OF-DFT and KS-DFT for N = 38 electrons in the main simulation cell.Additionally, we compare the results with the exact solution for the UEG with the LDA XC functional, where G LDA (k) ∼ k 2 is the local field correction of the UEG in the long wavelength limit.
To compute G LDA (k), we used the compressibility sum-rule [79], FIG. 4. Density profile along the perturbation direction at A = 0.5 and k ≃ 1.208k 0 F .The results are computed using OF-DFT with WT functional and KS-DFT and setting the XC functional to zero for both.The solid line is the density values n = n0 + δn with δn computed using the RPA density response function Eq. ( 14) of the UEG in Eq. ( 13).
For completeness, we note that a numerical analysis of Eq. ( 16) has been provided in the recent Ref. [56] both at ambient conditions and in the warm dense matter regime.
From Fig. 2, we observe that the WT based data is in perfect agreement with the exact solution (15).The same is the case for the KS-DFT data, which validates the solution (15).The MGP based results strongly deviate from the exact solution (15) and the KS-DFT data at 1 ≲ k/k 0 F ≲ 2. LMGP results follow the benchmark result despite showing small deviation at around 1.5 ≲ k/k 0 F ≲ 2. The PGSL (PGS) based OF-DFT data is in good agreement with the solution (15) and the KS-DFT results at k/k 0 F ≲ 2 (k/k 0 F ≲ 1).These relative deviations of the OF-DFT results from the exact data take place at exactly same wavenumbers as for the ideal electron gas case shown in Fig. 1a).Clearly, the inaccuracies in the OF-DFT results originate from the approximations in the KE functional.Therefore, we can isolate these errors from other possible effects by comparing the OF-DFT data and KS-DFT data computed with zero XC functional.Although it is trivial to see for the UEG, in a strongly inhomogeneous case, various types of numerical error cancellations might mask to a certain degree the deficiencies of a KE functional.

B. Strongly perturbed electron gas
Let us now consider a strongly perturbed electron gas.Specifically, we consider an electron gas with zero XC energy in order to analyze the errors due to the approximation of the KE functional in OF-DFT.To generate strong deviations from the uniform case, we use A = 0.1, A = 0.5, and A = 1.0.In Fig. 3, we show the values of the largest positive (δn ⊕ ) and negative (δn ⊖ ) deviations of the density from the mean value n 0 .We note that the |δn ⊖ | value is physically constrained to be smaller than n 0 since the electron density cannot have negative values.In Fig. 3, we present the results computed using the KS-DFT simulations with zero XC functional at wavenumbers of the external perturbation k min , 2k min , 3k min , and 4k min .The data is computed for N = 38 electrons.Correspondingly, we have k min = 2π/L ≃ 0.604k 0 F .The density perturbation values of less than ±0.1n 0 are indicated by the grey area.From Fig. 3, one can see that, at A = 0.1, we have |δn ⊕ | ≥ 0.1n 0 and |δn ⊖ | ≥ 0.1n 0 .These density deviations from the uniform case are already beyond the linear response regime [46,47] and can be characterized as a strongly inhomogeneous electron gas.A further increase in the amplitude of the external perturbation to A = 0. Due to the use of the local density approximation and various constrains in KE functionals (such as the positivity of the density and Pauli potential, and an exact single particle limit), OF-DFT is expected to describe strongly inhomogenous systems with a certain accuracy.This has motivated us to look at the density response of the electron gas to a strong external perturbation.In Fig. 5, we show the density profile computed for the perturbation amplitude A = 0.5 and the wavenumber k = 2k min us- ing the WT functional in OF-DFT and using KS-DFT with the XC functional set to zero in both.We compare the KS-DFT and OF-DFT results with the density values n = n 0 + δn from the linear response theory with δn being computed using the RPA density response function Eq. ( 14) of the UEG in Eq. (13).From this figure, we see that the linear response theory gives nonphysical negative values for the density since A = 0.5 is well beyond the weak perturbation regime [34,46].In contrast, the OF-DFT results based on the WT functional are in close agreement with the KS-DFT results and, as expected, always positive.We do not consider higher order response functions of the UEG appearing in the Taylor expansion of the density perturbation with respect to weak external perturbation [80,81].In fact, the considered KE functionals are designed to reduce (analytically at all or certain wavenumbers) to the Lindhard function in the UEG limit.Nevertheless, it can be seen from a corresponding Taylor expansion that the first order density response dominates over contributions of higher order terms.As we show next, this correlates with the observation that the KE functionals reproducing correct first order density response (at a given wavenumber of the perturbation) provide more accurate results.
The profile of the density distribution along the perturbation direction is shown in Fig. 4 for A = 0.1 and k = 2k min for both KS-DFT and OF-DFT simulations with different KE functionals.We observe that the OF-DFT data computed using the WT, LMGP, PGSL, and PGS KE functionals are in good agreement with the KS-DFT data.This clearly illustrates that the KE functionals based on the linear response function χ Lin (k) can remain valid beyond the weak density perturbation regime (given that they reproduce χ Lin (k) in the UEG limit).In contrast, the MGP and LKT based results show significant errors in the density values compared to the KS-DFT data.This is not surprising since these functionals are not able to accurately describe the ideal density response function of the UEG at the considered wavenumber k = 2k min ≃ 1.21k 0 F (see Fig. 1a)).To further analyze the quality of the considered KE functionals, we present histograms of the largest values of the error in the density accumulation and depletion regions at different values of the perturbation amplitude A and wavenumber k in Fig. 6.We consider perturbation amplitudes A = 0.1, A = 0.5 and A = 1 and wavenumbers k = k min ≃ 0.6k 0 F (panels a) and b)), k = 2k min (panels c) and d)), and k = 3k min (panels e) and f)).To quantify the error in the electron density, we measure the relative density deviation of the OF-DFT results relative to the KS-DFT data, We compute the ∆ n values for the density accumulation (δn(r) > 0) and density depletion (δn(r) < 0) regions separately.The former is presented in the panels labeled by ⊕ and the latter in the panels labeled by ⊖ in Fig. 6.The grey areas depict error values less than ±1%.From Fig. 6 we see that: The WT and LMGP KE functionals deliver rather accurate results with ∆ n < 5% for δn/n 0 ≤ 2.4 and all considered wavenumbers; The PGSL and PGS KE functionals are accurate at k ≲ k 0 F for the harmonically perturbed inhomogeneous electron gas with δn/n 0 ≲ 0.5; The LKT KE functional is accurate at k ≲ 0.5k 0 F for the harmonically perturbed inhomogeneous electron gas with δn/n 0 ≲ 0.5; The MGP KE functional is accurate at k ≲ 0.5k 0 F for the harmonically perturbed inhomogeneous electron gas with δn/n 0 < 0.5 and cannot be considered as reliable at larger wavenumbers.
Summarizing the key findings from the data presented in Fig. 6 and for the UEG limit, we conclude that the accuracy of a KE functional for the inhomogeneous regime strongly correlates with its accuracy for the density response in the limit of the UEG.

C. Bulk Aluminum with fcc lattice
As an example for the application of our approach to a real system, we next consider bulk Al with an fcc structure.In general, the electronic density in materials is inhomogeneous.This is illustrated in Fig. 7a), where we show the electron density in the unperturbed case using volume rendering.The results presented in Fig. 7 are computed using KS-DFT with an LDA [82] XC functional.Additionally, we show the density perturbations created by an external harmonic field with the amplitude A = 0.01 and the wavenumbers b) k ≃ 0.9k 0 F , c) k ≃ 1.8k 0 F , and d) k ≃ 2.7k 0 F .In contrast to the UEG case, one can see from Fig. 7 that the density perturbation is also inhomogeneous in the direction perpendicular to the external harmonic perturbation wave vector (which is indicated in Fig. 7 by an arrow).In order to quantitatively analyze the difference between OF-DFT results and KS-DFT results, we consider the projection of the density and density perturbation onto the x-axis along which the harmonic perturbation is directed.
In Fig. 8, we present the corresponding results for the density projection along the x-axis computed using OF-DFT with different KE functionals and using KS-DFT.We used the same LDA XC functional in all calculations.The density value is shown in the units of the mean density of the valence electrons n 0 ≃ 1.81 × 10 23 cm −3 .
Examining the data in Fig. 8, we see that a meaningful analysis of the OF-DFT results can be performed by looking at ∆n = n − n 0 , which is the deviation of the density profile from the mean electron density n 0 ; the latter is depicted as the dashed horizontal line in Fig. 8. First, we note that, in the case of the KS-DFT data, ∆n KS−DFT has the largest value about 0.075 n 0 , i.e., about 7.5% of the mean density.As expected for the valence electrons of metals, this indicates that the valence electrons are weakly perturbed by the ions in the bulk region of Al with an fcc lattice.We note that the largest deviation amplitude in the density accumulation region is smaller than the largest density deviation magnitude in the density depletion region by about 17%.This means that the density inhomogeneity is not symmetric with respect to the positive and negative deviations from the mean density.Comparing the ∆n data computed from the OF-DFT simulations with the results from the KS-DFT simulations, we find that the WT and LMGP based results closely reproduce the KS-DFT data, FIG. 9.
Projection of the density perturbation profile δn(x)/n0 on the x-axis in units of the mean valence electron density (n0 ≃ 1.81 × 10 23 cm −3 ).The results are computed for the primitive unit cell of fcc Al using OF-DFT with different KE functionals and KS-DFT; an LDA XC functional has been used for all calculations.The density perturbation is induced by an external harmonic field with the amplitude A = 0.01 and wavenumbers a) with a maximum relative error of about 3%.The MGP data exhibit a larger disagreement with the KS-DFT data with a maximum relative error in ∆n of about 9%.In the case of the PGSL KE functional, the largest relative errors in ∆n are about 35%.For the LKT and the PGS KE functionals, we observe the largest relative errors in ∆n about 44% and 47%, respectively.Therefore, as one might expect, the semi-local GGA KE functionals are less accurate for the OF-DFT calculations of the density compared to other considered fully nonlocal KE functionals in this case.The disagreement between OF-DFT and KS-DFT results are most pronounced around the maxima and minima of the density.We do not consider the effect of these observations on the calculation of other bulk properties since the effect of density errors can be masked or canceled to a certain degree by other effects, e.g., the behavior of the total energy [39].We next proceed to the analysis of the results for the density response to the external harmonic perturbation computed using OF-DFT and KS-DFT.In Fig. 9, we show the density perturbation δn = n A − n A=0 induced by an external harmonic potential with the amplitude A = 0.01 and wavenumbers a) and c) k ≃ 2.7k 0 F .Additionally to the DFT data, we provide the density perturbation values computed using the density response function of the UEG according to Eq. ( 15) at n 0 ≃ 1.81 × 10 23 cm −3 (dashed line).Comparing the KS-DFT results with the UEG model results, we see that the UEG model is able to reproduce the KS-DFT data at all considered wavenumbers despite of the microscopic inhomogeneities in both the unperturbed and perturbed density distributions (as illustrated in Fig. 7).To understand this observation, we note that the used harmonic perturbation can be considered as macroscopic since it acts at all space points and the reaction of the system to this perturbation observed in Fig. 9 is a macroscopic density response since we have performed averaging over the density values in the direction perpendicular to the perturbation wavevector.The thus obtained macroscopic static density response of the valence electrons in metallic Al is accurately described by the UEG model.Although it is considered to be common knowledge, as far as we know, this had not been demonstrated quantitatively in this way.This example demonstrates also how the direct perturbation approach can be used to acquire a physical picture of electronic properties in real materials.
Let us now summarize the performance of the considered KE functionals within OF-DFT for the description of the density perturbations presented in Fig. 9. Despite being generated on top of an inhomogeneous density, the quality of the density response from OF-DFT relative to the KS-DFT results correlates with the observations for the UEG discussed in Sec.IV A. The WT and LMGP based results closely reproduce the KS-DFT data.In spite of the inaccurate unperturbed density values, the PGSL based results for the density perturbation by the weak external harmonic field are in good agreement with the KS-DFT data at k = k min ≃ 0.9k 0 F , k = 2k min ≃ 1.8k 0 F , and k = 3k min ≃ 2.7k 0 F (being slightly worse for the latter).This can be due to the aforementioned effect of the averaging of the density inhomogeneities and due to the small magnitudes of these density inhomogeneities compared to the mean density value.The PGS data is similar to that of PGSL at k = k min and k = 2k min , but is not able to adequately describe the density perturbation at k = 3k min .In the case of the MGP functional, we see from Fig. 9 that the corresponding results have large deviations from the KS-DFT data at wavenumbers k ≃ 0.9k 0 F and k = 3 ≃ 2.7k 0 F , and are in close agreement with the KS-DFT data at k ≃ 1.8k 0 F .The LKT based results significantly underestimate the density perturbation values at k = 2k min and are in rather close agreement with the KS-DFT data at k = k min and k = 3k min .This is consistent with the pattern we have observed considering the density response of the UEG in Sec.IV A.

D. Bulk Silicon with fcc, bcc, and cd lattice
Next, we consider bulk Si with an fcc lattice, bodycentered cubic (bcc) lattice, and semiconducting crystal diamond (cd) phase.In Fig. 10 and Fig. 11, we show the density perturbation values at k = k min = 2π/L, k = 2k min , and k = 3k min , where L is the corresponding length of the unit cell.We set the XC functional to LDA [82] in both KS-DFT and OF-DFT simulations.For the fcc lattice, we have k min ≃ 0.8k 0 F and for the bcc lattice we have k min ≃ 1k 0 F .Similar to the fcc Al case and for the same reasons, we see that the density perturbations in the fcc and bcc Si are accurately described by the UEG model.Correspondingly, the WT functional based OF-DFT data is in close agreement with the KS-DFT data for both fcc and bcc structures.Similarly to the case of the fcc Al, the performance of the other considered KE functionals strongly correlates with their accuracy for the density response in the limit of the UEG at relevant wavenumbers.
In the case of the cd Si phase, we see from Fig. 12 that the UEG model is less accurate.This is expected since cd Si is a semiconductor.The disagreements between the UEG model and KS-DFT results are particularly pronounced at k = 2k min and k = 3k min (for cd Si we have k min ≃ 0.64k 0 F ).Despite being based on the UEG model (in the sense of Eq. ( 1)), the WT, LMGP, PGSL, and PGS KE functionals (used in the OF-DFT simulations) provide better agreement with the KS-DFT data than the UEG model.This demonstrates the versatility of these functionals for the description of bulk metals and semiconductors [13,[39][40][41][42][43].Additionally, we observe that the LKT and PGS KE functionals provide results similar to each other for cd Si.

V. CONCLUSIONS AND OUTLOOK
We have demonstrated the application of the direct perturbation approach for the test and analysis of various KE functionals for OF-DFT calculations.By measuring the density response of the UEG to an external harmonic perturbation, we are able to crosscheck whether a KE functional satisfies the constraint (1).We emphasise that this tool is independent from the analytical reasoning used to construct a given KE functional.As a demonstration, we analyzed the density response generated by OF-DFT simulations using the nonlocal WT, MGP, LWT, LMGP functionals, a Laplacian-meta-GGA level PGSL functional, and GGA level LKT and PGS functionals.The WT, MGP, LWT, LMGP functionals were built utilizing the constraint (1) for all wavenumbers and the PGSL (PGS) is designed to respect the constraint (1) at k < 2k 0 F (k < k 0 F ).However, using the direct perturbation approach, we found that the MGP KE functional violates the constraint (1) at intermediate wavenumbers 0.2 ≲ k/k 0 F ≲ 2.5.This illustrates the utility of the direct perturbation approach for testing KE functionals.
Going beyond the UEG limit, we analyzed the results computed using the considered KE functionals for the harmonically perturbed inhomogeneous electron gas over a wide range of wavenumbers and density perturbation degrees.We found a strong correlation between the performance of the KE functionals in the UEG limit and in the strongly inhomogeneous case.This empirically demonstrates the importance of the constraint (1) for the construction of accurate KE functionals.Furthermore, for the example of the PGSL and PGS KE functionals, we numerically validated the mapping of s and q-constructed using density gradients-on k/(2k 0 F ) and k 2 /(2k 0 F ) 2 in the long wavelength expansion of the inverse density response function for the construction of KE functionals.This is an important result since such a mapping is also used for the construction of XC functionals using the local field correction (XC kernel) of the UEG (e.g., see Ref. [83]).
The application of the direct perturbation approach for the analysis of the density response properties and performance of the KE functionals for real materials are demonstrated for the example of Al with an fcc lattice structure and Si with fcc, bcc, and cd phases.We demonstrated that the comparison of KS-DFT data for the density perturbation to the UEG model allows one to understand the role of the microscopic density inhomogeneity (induced by ions) with respect to the density response properties of the valence electrons.Despite of the microscopic inhomogeneities present in Al and in Si with fcc and bcc structures, the macroscopic density response of the valence electrons is accurately described by the UEG model.As the result, the quality of the considered KE functionals for the description of the density perturbation in fcc Al, fcc Si, and in bcc Si is similar to that for UEG case.Despite being based on the UEG model, the WT and LMGP nonlocal functionals provide an adequate description of the density response of cd Si state.For this case, the PGS and LKT have similar quality and the PGSL performs best among considered semi-local functionals.
In summary, the results clearly show that the quality of the static response functions of periodic bulk systems (showed for the example of the UEG, fcc Al, fcc Si, and bcc Si) require KE functionals that specifically encode UEG response behavior in the limiting case of constant ground state densities.While such a requirement is not important for perturbations with small wavevectors, it is crucial for high-k perturbations with repercussions beyond ground-state OF-DFT.
When modeling optical and ultrafast electronic properties of materials, time-dependent OF-DFT simulations of bulk systems, therefore, will need to move away from semilocal KE functionals and instead employ fully nonlocal functionals to be able to provide physical results for high-k perturbations for periodic solids.
For modeling warm dense matter, Graziani et al [20] have recently pointed out the possible importance of quantum nonlocality effects for the description of the shock propagation.Preliminary results computed using a vW functional based quantum KE potential indicate that the induced density change at a shock front can reach about 0.5n 0 for wavenumbers k ≲ k 0 F .Taking into account these findings and the result of the present analysis, we suggest that the potential generated by a semilocal (e.g., PGSL and PGS) KE functional can be used for a more reliable investigation of the impact of quantum nonlocality on the shock propagation in warm dense matter with strongly degenerate electrons.We single PGSL and PGS KE functionals due to their relatively simple form compared to the fully nonlocal KE functionals.This is advantageous for the implementation into existing hydrodynamics codes for which it is important to minimise a computational overhead due to calculation of the force field.Additionally, semilocal KE functionals such as PGSL are able to provide correct density response properties of electrons at k ≲ k 0 F , which is important for the adequate description of the ion screening in warm dense matter [54,84,85].We note that the finite temperature LKT functional [86] and finite temperature WT functional [53] can be used to extend the present analysis into the WDM regime with partially degenerate electrons, and to the simulation of shock propagation at high temperatures relevant for inertial confinement fusion experiments.
In this work, we considered the KE functionals constructed using the Lindhard density response function.There are other types of the density response function that can be used for the construction of a KE functional, such as the quadratic density response function of the ideal electrons gas [87,88].Since the analytical forms of these density response functions are known, the presented harmonic perturbation based approach can be used as an independent tool for testing the correctness of the implementation of the KE functionals built using any of these density response functions.
Finally, we also note that a more detailed study of the response of semiconductors is the subject of future work, where in addition to the nonlocal KE functionals with density-independent kernel considered here, one can employ KE functionals based on the jellium-with-gap model [45].

FIG. 1 .
FIG. 1. Linear density response function of the UEG in the ground state at metallic density a) with screening (RPA) and b) without screening.The results are computed using OF-DFT with different KE functionals and KS-DFT setting the XC potential to zero.The solid line represents the exact analytical results.The simulation data are computed for N = 7168, N = 66, N = 38 and N = 14 electrons in the simulation cell.

FIG. 2 .
FIG. 2. Density response function of the correlated UEG in the ground state at metallic density computed using OF-DFT with different KE functionals and KS-DFT with the LDA XC functional.The dashed line represents the analytical result Eq. (15).The results are computed for N = 7168 electrons in the simulation cell.

FIG. 5 .
FIG. 5. Density perturbation profile δn(r) along the perturbation direction at A = 0.1 and k ≃ 1.208k 0 F .The results are computed using OF-DFT with different KE functionals and KS-DFT setting the XC potential to zero.

FIG. 6 .
FIG. 6.The largest value of the error in the relative density deviation of the OF-DFT data from the reference KS-DFT data, Eq. (17), in the region with δn(r) > 0 (top panel, denoted as ⊕) and in the region with δn(r) < 0 (bottom panel, denoted as ⊖) for A = 0.1, A = 0.5, and A = 1 at k ≃ 0.604k 0 F .The grey area indicates |∆ n⊕| ≤ ±1 % and |∆ n⊖| ≤ ±1 %.The results are computed for the electron gas with XC potential set to zero.

FIG. 7 .
FIG. 7. Bulk electron density of Al with an fcc structure.Volume rendering is used to visualize a) the electron density for the unperturbed case (A = 0), and the density perturbation due to an external harmonic perturbation with the amplitude A = 0.01 and wavenumbers b) k ≃ 0.9k 0 F , c) k ≃ 1.8k 0 F , and d) k ≃ 2.7k 0 F .Al atoms are shown as green spheres.

FIG. 8 .
FIG. 8.Profile of the density projection n(x)/n0 on the x-axis in units of the mean valence electron density (n0 ≃ 1.81×10 23 cm −3 ).The results are computed for the primitive unit cell of fcc Al using OF-DFT with different KE functionals and KS-DFT; an LDA XC functional has been used for all calculations.

FIG. 10 .FIG. 11 .
FIG. 10.Density response δn(x)/n0 in units of the mean valence electron density (n0 ≃ 2.82 × 10 23 cm −3 ).The results are computed for the primitive unit cell of fcc Si using OF-DFT with different KE functionals and KS-DFT; an LDA XC functional has been used for all calculations.The density perturbation is induced by an external harmonic field with the amplitude A = 0.01 and wavenumbers a) k ≃ 0.8k 0 F , b) k ≃ 1.61k 0 F , and c) k ≃ 2.42k 0 F .

FIG. 12 .
FIG. 12. Density response δn(x)/n0 in units of the mean valence electron density (n0 ≃ 2.03 × 10 23 cm −3 ).The results are computed for the primitive unit cell of cd Si using OF-DFT with different KE functionals and KS-DFT; an LDA XC functional has been used for all calculations.The density perturbation is induced by an external harmonic field with the amplitude A = 0.01 and wavenumbers a) k ≃ 0.64k 0 F , b) k ≃ 1.28k 0 F , and c) k ≃ 1.92k 0 F .

FIG. 13 .
FIG. 13.Density profile along the perturbation direction at A = 0.01 and k ≃ 1.208k 0 F .The OF-DFT results are computed using the WT functional.The solid line represents the functional form χ(k)2A cos(kx) with χ(k) being computed from a fit to the OF-DFT results.