Nonlocal kinetic energy functionals by functional integration Free

Since the seminal work of Kohn and Sham,¹ the energy functional of the electron density can be written in terms of the noninteracting kinetic energy (KE) functional, T_s[ρ], and the associated exchange-correlation (xc) functional, E_xc[ρ], as follows:

The nuclear-nuclear interaction, E_NN, is density independent. The electron-nuclear interaction, E_eN[ρ] = ∫$vext$(r)ρ(r)dr, is linear in the density, and the other terms [for closed-shell systems and introducing the Kohn–Sham (KS) orbitals, {ϕ_i[ρ](r)}, including the yet unknown exchange–correlation (xc) energy density, ε_xc[ρ](r)] take the form

There are two different and theoretically equivalent ways to find the electron density associated with a given external potential. One prescribes the direct minimization of the density-functional theory (DFT) Lagrangian,

which is known as Orbital-Free DFT (OF-DFT). The other way (according to Kohn and Sham¹) prescribes solving the following Schrödinger equation:

If one decides to carry out the density search according to the minimization of the DFT Lagrangian in Eq. (5), then there is no need to invoke the concept of Kohn–Sham (KS) orbitals. Instead, the electron density, ρ(r), can be utilized as the only variational parameter. The drawback is that in addition to approximations needed for the xc functional, in OF-DFT, approximations for the kinetic energy density functional (KEDF), T_s[ρ], are also needed.

Local and semilocal (i.e., dependent on the value of the electron density and its gradient at a point in space) parametrizations of KEDFs have shown potential^2–4 and displayed a more than satisfactory agreement with KS-DFT for simulations of warm dense matter.⁵ Recent work by the Della Sala group⁶ advocated for employing the Laplacian in addition to the gradient in the formulation of T_s[ρ] and the possibility of reproducing accurate kinetic energy densities.⁷ 

Nonlocal versions of T_s[ρ], such as the Wang-Govind-Carter (WGC),⁸ the one proposed by Perrot,⁹ the latest Huang-Carter (HC),¹⁰ and others,¹¹ typically improve over local and semilocal functionals,^10,12–16 and all share the following form:

where $TTFρ$ is the Thomas-Fermi kinetic energy,^17,18$TvWρ$ is the gradient correction according to von Weizsäcker,¹⁹ and $TNLρ$ is the remaining contribution which is by construction of nonlocal character. The general form of T_NL is

where the function $ωTNL$ is commonly called “kernel of the nonlocal KEDF” and α and β are the suitable exponents.

Unfortunately, the current state of the art for OF-DFT is that, although the algorithms have dramatically improved over the years and are typically orders of magnitude faster than KS-DFT, their application to semiconductors and molecular systems is still problematic.^10,12,13,15 For example, WGC^8,15,20 and Wang-Teter (WT)²¹ KEDFs can reproduce remarkably accurate bulk properties for light metallic systems, but are less accurate for semiconductors and nonsimple metals. The HC KEDF¹⁰ produces generally reliable bulk properties for semiconductors. However, it underestimates the electron density in the bonding regions of crystal diamond (CD) silicon and III-V semiconductors, it is less accurate than WT for metals, and it is substantially more computationally expensive than any other KEDFs.^10,22 Density decomposition schemes that combine nonlocal with other functionals have found applicability to molecules, materials, and alloys^12,22,23 and thus appear to be an interesting avenue of research.

In this work, we aim at formulating a KEDF that is (1) at least as accurate as WGC and WT KEDFs for nonsimple metals such as the body-centered cubic (BCC) and face-centered cubic (FCC) phases of silicon; (2) at least as accurate as the HC KEDF in predicting bulk properties of CD silicon and III-V semiconductors, especially in regard to the self-consistent electron density in bonding regions; and (3) computationally inexpensive.

It is known that the exact T_s functional satisfies a number of mathematical relations also known as exact conditions. It is arguable²⁴ that the more exact conditions the functional satisfies, the closer the functional necessarily is to the exact functional. It is through the imposition of exact conditions that the Trickey group developed generalized gradient approximation (GGA) KE functionals approaching the accuracy of KS-DFT in the limit of high temperature.⁵ By imposing exact conditions on the asymptotic behavior of the linear response function, the Carter group formulated the HC KEDF which could be applied to systems with record inhomogenous density for OF-DFT simulations.¹⁰ 

Let us enumerate three exact conditions that we believe are particularly important: (1) existence of the “Kinetic electron,”²⁵ (2) hypercorrelation, i.e., potentials are related by functional integration to their respective second functional derivatives, and (3) the Lindhard response in the Free Electron Gas (FEG) limit.

The first unconstrained functional derivative of the DFT Lagrangian with respect to the density yields the Euler equation of DFT.²⁶ Namely,

By taking the negative of the Laplacian of each term and rearranging, $−∇2vTs[ρ](r)=∇2vH[ρ](r)+∇2vxc[ρ](r)+∇2veN(r)$⁠. After solving the associated Poisson equations for each term in the long-range limit,²⁷ we find

where ρ_N is the charge density of the nuclei (ions) and ρ_xc is related to the exchange-correlation hole. At long ranges, ρ_xc is identical to the xc hole.²⁷ Thus, at long ranges, $ρTs$ is the kinetic energy equivalent of the xc hole. However, upon simplification of the 4π terms and integration over all space, we find that

where the integration of the negative of the xc hole is 1, the integration of the electron density is N_e, and the integration of the nuclear density is N_N, and the total charge of the system is defined as N_N − N_e.

We define $ρTs$ as the Kinetic electron, as for neutral systems, it integrates to +1 in contrast to the xc hole which integrates to −1.

It is possible to recover an entire functional from the corresponding functional derivative using the following line integral:^28,29

where we have chosen the linear interpolation path which we describe below.

As the potential is the functional derivative $v(r)=δE[ρ]δρ(r)$⁠, the result of the integration is path-independent. The scaled density, ρ_t, can then have almost any form,²⁸ the simplest of which is given by a linear interpolation between the vacuum and the density,

Thus, the t-integral equation reduces to

The previous integral defines the energy density in terms of the hypercorrelated potential,

Hypercorrelation is a term coined by Burke³⁰ to emphasize the fact that the energy density is directly related to the physical (correlated) potential evaluated for an array of electron densities ranging from the vacuum to the true density.

The above equation can be considered yet another exact constraint any functional (including the kinetic energy) should satisfy.

Hypercorrelation does not only concern the energy density but also the potential. That is, the potential must be related to the second functional derivative by

It is important to notice that in the definition of the scaled second functional derivative on the rhs, ρ_t appears only on one of the two functional derivatives.

In light of Eq. (16), Eq. (14) can be recast in terms of the second functional derivative as

We should remark that due to the derivative discontinuity,³¹ the line integration in the above equations will cross discontinuities in the integrands when integrating the kernel to yield a potential and also when integrating the potential to obtain an energy value. These discontinuities are integrable and do not pose formal problems. Equation (17) may worry the reader because, in principle, terms proportional to the Dirac delta function in real space could appear in the second functional derivative and could lead to density linear terms in the potential and density quadratic terms in the functional. This of course is not possible in theory due to the known density scaling relations of the exact KEDF functional.³² In practice, this also does not occur because the chosen second functional derivative (vide infra) is such that these terms do not appear. Similar arguments lead to the realization that density linear terms should not appear in the KEDF.

An exact condition for a KE functional is its relationship to the linear response function of the FEG. There is an established relation between the KE functional and the linear response function.^33–35 Taking a functional derivative of each term in the Euler equation (9),

and the response function is given by $χs(r,r′)=δvs(r)δρ(r′)−1$⁠. Unfortunately, the above relationship is only valid at self-consistency, and trying to impose it directly would lead to impractical algorithms. However, for the FEG, we have a simplified relationship, χ_s(r, r′) = χ_s(|r − r′|), resulting in a one-dimensional Fourier transform, χ_Lind(η), the Lindhard function,^36,37

with $η=q2kF$⁠, and the Fermi wavevector is given by $kF=(3π2ρ0)1/3$⁠, with $ρ0=NeV$⁠.

Thus, in the limit of a constant and periodic electron density, we must impose

where $F^⋅$ stands for the Fourier transform in q, the conjugate variable of |r − r′|. This exact condition is at the core of the original Hohenberg and Kohn functional³⁵ which inspired the entire class of nonlocal KE functionals, well summarized in Ref. 33. The function G_Lind(η) is introduced here as

We cast our developments in terms of Eq. (8). Particularly, in the following, we will only focus on the nonlocal term of the KEDF. We begin by espousing the idea that in the limit of modeling the Free Electron Gas (FEG), a nonlocal KEDF should recover the exact FEG linear response function. This can be imposed by hypercorrelation, Eq. (16). Namely,

where $δvTNL[ρt](r)δρt(r′)=δ2TNL[ρt]δρ(r)δρt(r′)$ will be dependent on ρ_t and consequently on t. Applying a local density approximation (LDA) on the polynomial terms of Eq. (20), $δ2TNLδρ(r)δρt(r′)$ can be rewritten as

where, following the prescriptions of Eqs. (7), (19), and (20), the Fourier transform of G_NL(|r − r′|) is taken to be

In the above, we have substituted ρ₀ → tρ₀, which gives the following Fermi wavevector $kF(t)=t133π2ρ013$⁠, and the reciprocal space variable is $η(q,t)=q2kF(t)=q2t133π2ρ013$⁠.

Equation (23) is derived from Eq. (22) in the following two steps: First, the G_NL term depending on the reciprocal space variable η is represented in terms of its conjugate real-space variable |r − r′|. In a second step, the Fermi wavevector multiplying G_NL has a $ρ0−1/3$ dependence that can be split into $ρ0−1/6⋅ρ0−1/6$ and, after symmetric LDA approximation, leads to the ρ(r)^−1/6·ρ(r′)^−1/6 dependence of the second functional derivative.

Thus, substituting Eq. (23) into Eq. (22), the nonlocal kinetic potential, $vTNL[ρ](r)$⁠, can be rewritten as

Once again, we point out that the integral above should cross discontinuities in the integrand as the density scaling parameter, t, is varied. This does not directly affect the integral above because the Lindhard function is continuous across all values of t.

Finally, the nonlocal kinetic potential can be obtained by

where the kernel in reciprocal ω_T(q) is obtained from integration by parts,

The new kernel, $ωTNL$⁠, is given by the kernel of the WT KEDF²¹ plus a correction term. We should point out that this correction term derived from assuming that the inverse Lindhard function is dependent on the electron density (this allowed us to include in the line integral the average density ρ₀). Thus, we expect $ωTNL$ in Eq. (27) to be the most appropriate for substantially nonuniform densities. For light (simple) metallic systems, the distribution of electron density is very close to the uniform electron gas, and thus the WT kernel is already very good for these systems. We will see in our pilot calculations that for more complex systems, such as (nonsimple) metallic phases of silicon, IV, and III-V semiconductors, the correction term in conjunction with an additional correction arising from the presence of the “Kinetic electron” remarkably improves upon the performance of WT.

We also note that in going from Eqs. (22) to (23), we have assumed that the inverse Lindhard function is a function of t multiplied by ρ₀. It is arguable that symmetry properties of the second functional derivative should be imposed. In the supplementary material, we carry out a detailed analysis of symmetrization of the kernel. There, we propose an arithmetic symmetrization as well as a geometric symmetrization of the kernel. We note that in either cases, the resulting kernels yield results that are equivalent to the ones obtained with the kernel in Eq. (27) (vide infra).

Equation (10) implies that the KE potential in reciprocal space is given by

where the tilde indicates the quantities which have been Fourier transformed in reciprocal space (the Fourier transform is kept one dimensional for the sake of simplicity). Unfortunately, the potential derived by Eq. (26) is zero in the low q limit (or at long ranges in real space). That is,

and with that the Kinetic electron is absent.

Thus, to impose a nonzero Kinetic electron, we decided to model it with a simple Gaussian multiplied by the square of the error function centered at q = 0 in order to remove the Coulomb singularity. This results in a new kernel, ω_MGP, obtained simply by modifying the expression for the kernel in Eq. (27). Namely,

Although modeling the Kinetic electron by a simple Gaussian function is a gross approximation, in the future, we commit to explore more accurate expressions inspired by the established work of others.^38,39 We are aware that the exchange hole (and thus also the Kinetic electron) is generally not spherical,⁴⁰ and thus in the future, we will also consider nonspherical parametrizations. In Fig. 1, four different kernels are plotted. It is clear that the presence of the Kinetic electron imposed through Eq. (30) affects dramatically the low-q limit of the kernel. The kernel resulting from our manipulations is smooth, which will result in smooth nonlocal kinetic energy potentials.

To understand the effect of the Kinetic electron, let us consider a very narrow Gaussian (b → +∞). This results in a spatially extended real-space Kinetic electron. This is the preferred shape for materials with a finite gap, in which the dielectric screening is small. The opposite case (b → 0) is preferred for metals, as the dielectric screening is large and the real-space spatial extension of the Kinetic electron is reduced. Thus, we expect that small b values are suitable to model metallic systems and large b values are suitable to model semiconductors. These observations are consistent with the finding^41–44 (also exploited by Huang and Carter¹⁰) that, for q → 0, the static response function for systems with a gap, χ → −Constant × q². An in-depth analysis of the relation between the low-q behavior of χ and the long-range properties of the potential will be the subject of a follow-up study.

We implemented the new functional resulting from the kernel in Eq. (30) according to the framework outlined in Eq. (26). The Mi-Genova-Pavanello (MGP) potential and energy functional are defined as

where $F^−1$ indicates the inverse Fourier transform and $vTF$(r) + $vvW$(r) is the local part of the kinetic energy potential composed of the sum of Thomas–Fermi and von Weizsäcker terms. To avoid spurious dependencies from the parameter a, we choose such boundary conditions as ω_MGP(q = 0) = 0.

From Eqs. (32) and (17), the energy expression is derived. Namely,

The simplicity of the above equation results from double integration of the kernel in Eq. (23). A similar expression would be recovered for the arithmetically symmetrized kernel, while a more complex expression involving a double integration of the inverse Lindhard function would be required for the geometrically symmetrized kernel.

We mention here that the so-called LQ and HQ functionals¹¹ (their potentials are constructed ad hoc to exactly reproduce the high-q or the low-q limits) feature a nonlocal kinetic energy that is also evaluated with a line integral similar to Eq. (32). The line integral in Ref. 11, however, is carried out only for the nonhomogeneous component of the density [i.e., ρ(r) − ρ₀].

The OF-DFT calculations are carried out with both modified versions of Ab initio orbiTaL-free density functionAl theory Software (ATLAS)^45,46 and PRinceton Orbital-Free Electronic Structure Software (PROFESS) 3.0.⁴⁷ Calculations involving the HC KEDF are performed with PROFESS 3.0, and all others are carried out in ATLAS and double checked with PROFESS 3.0 for consistency. The benchmark KS-DFT calculations are carried out with the same local pseudopotentials as the OF-DFT simulations. We employed the optimal effective local pseudopotentials (OEPP)⁴⁸ and the bulk-derived local pseudopotentials (BLPS) of Huang and Carter.⁴⁹ KS-DFT calculations with OEPP are performed with Quantum-ESPRESSO (QE),⁵⁰ while the ones employing BLPS are carried out with ABINIT.⁵¹ In all calculations, the LDA xc functional by Perdew and Zunger⁵² is adopted.

Crystal diamond (CD), body-centered cubic (BCC), and face-centered cubic (FCC) silicon phases, as well as nine III-V cubic zincblende (ZB) semiconductors, are selected as benchmark systems. We compute total energy curves as a function of the lattice constant and extract the cell volume, V₀, the minimum energy, E₀, and the bulk modulus, B₀. These were computed using the prescription of Huang and Carter,¹⁰ fitting the energy curves vs volume against Murnaghan’s equation of state⁵³ within 0.95 V₀ to 1.05 V₀.

In all OF-DFT calculations, the grid spacing of 0.2 Å in ATLAS and kinetic energy cutoffs 1600 eV in PROFESS 3.0 are sufficient for well-converged total energies (1 meV/cell). In KS-DFT simulations, for bulk properties, the 20 × 20 × 20 k-point meshes and 800 eV kinetic energy cutoffs are used. To obtain smooth electron density, a denser grid (larger kinetic energy cutoffs) is adopted for both OF-DFT and KS-DFT calculations to keep the real space electron density as represented on 54 × 54 × 54 for ZB GaAs and 36 × 36 × 36 for CD silicon.

We have benchmarked MGP against KS-DFT and compared it to WT, WGC, and HC KEDFs for metallic and semiconducting phases of Si, as well as the ZB phase of nine II-V semiconductors. In the following, we will only compare the so-called “optimal” parameters for all KEDFs, including MGP. We should remark that the optimal parameters are pseudopotential dependent. Thus, we optimize all KEDF parameters accordingly. With exception of WT, all other KEDFs have two adjustable parameters: MGP has a and b, defining the Kinetic electron term; WGC has the effective electron density, ρ^*, and the averaged Fermi wavevector parameter, γ; and HC has β and λ, defining the long-range portion of the functional.

The optimal parameters for MGP functionals along with the optimal γ and ρ^* for WGC, optimal β and λ for HC functionals, and calculated quantities for all KEDFs are collected in Table I.

Since BCC and FCC structures of Si are metallic, we expect that WGC KEDF produces more or comparably accurate bulk properties compared to HC KEDFs. Thus, for these metallic phases, we just compare MGP results to WGC, WT, and the benchmark KS-DFT results. Comparing to KS-DFT, MGP reproduces bulk properties overall improving on WT and WGC KEDFs with both OEPP and BLPS. This is especially the case for the total MGP energies which deviate from KS-DFT by less than 5 meV/cell in all cases.

Modeling CD Si with OF-DFT has historically been a challenge which has been addressed by several nonlocal functionals with a density-dependent kernel.^8,10,20,22 For example,¹⁵ WT is unable to reproduce a bound curve for CD Si. As MGP functional has a density-independent kernel, we had no expectations that CD Si would be modeled correctly. Table I, however, shows that MGP is capable of producing bound energy curves and overall bulk properties that are close to the KS-DFT benchmark compared to other KEDFs. To the best of our knowledge, MGP is the only KEDF with density-independent kernel capable of simultaneously reproducing KS-DFT equilibrium total energies, bulk moduli, and equilibrium volumes for both metallic and semiconductor phases of silicon.

We also tested a modified version of MGP only including WT and the Kinetic electron term. We found that the modified functional could not reproduce the benchmark despite efforts in optimizing the a and b parameters. To give an example, we have tried to optimize parameters a and b for the CD phase of silicon and found the best values to yield a reasonable V₀, E₀ = −217.427 eV, and B₀ = 126 GPa. These values are severely deviated from the KS-DFT benchmark. Thus further analysis of the WT+KE model is unnecessary. This indicates that, to achieve accurate results, it is not enough to simply add the correction terms separately. Instead, the correction term derived from functional integration and the Kinetic electron must be included together for MGP to approach both metallic and semiconducting phases.

In this section, we present the benchmark MGP against KS-DFT for binary III-V semiconductors in the ZB phase. These are challenging systems for KEDFs. For example, it was shown^10,12,15 that WGC is not appropriate and that the proper long-range behavior of the kernel needs to be included. We have developed MGP with semiconductors and finite systems in mind, and specifically, the Kinetic electron term aims at correcting the KEDF kernel in its long-range (low q) behavior.

Table II lists the performance of MGP compared to HC as well as KS-DFT. The results indicate that both MGP and HC functionals can reproduce accurate equilibrium volume and bulk moduli. The predicted equilibrium volumes with MGP and HC lie within 2% of the KS-DFT results for all considered semiconductors. For the bulk moduli, MGP and HC are within 10 GPa compared to KS-DFT. MGP total equilibrium energies are within 5 meV/cell to the KS-DFT benchmark. This improves upon HC which finds itself within 42 meV/cell. Overall, however, the quality of the calculated MGP energies is similar to HC for the considered systems.

We have also evaluated MGP by comparing the ground state electron densities produced by OF-DFT with MGP and HC KEDFs, as well as KS-DFT for CD Si and ZB GaAs. Once again, we employ optimal parameters for the energy functionals. Figures 2 and 3 display the electron density distribution along the [111] direction for CD silicon and ZB GaAs, respectively. Inspection of the figures reveals that the MGP electron density reproduces KS-DFT quite well. While in both cases HC tends to underestimate the density in the bonding region, MGP is much closer to KS-DFT. Conversely, in the space between non-covalently bonded atoms, HC slightly underestimates the KS-DFT density, while MGP slightly overestimates it.

We have tested several functionals for their performance in reproducing KS-DFT electron densities for these systems and found that WT does a remarkably good job. In this respect, because MGP is based on WT, it inherits some of its good traits.

To evaluate the computational efficiency of MGP, we compare the computational cost for optimizing the electron density in CD Si (2 atoms), 2 × 2 × 2 (16 atoms), 4 × 4 × 4 (128 atoms), and 8 × 8 × 8 (1024 atoms) supercells with WGC, WT, and MGP KEDFs, respectively. All calculations are performed on a single thread with PROFESS employing the BLPS using a large kinetic energy cutoff (4000 eV). For all KEDFs, we employ optimal parameters (WGC: γ = 4.2, ρ^* = ρ₀; WT: α = β = 5/6; MGP: a = 0.364 and b = 0.570). Since the computational cost of HC KEDF is much higher than that of other functionals (about 1000 times higher than MGP for CD Si),^10,22 we exclude it in the comparison. In principle, the only difference in the computational cost between MGP and WT KEDFs is the initial kernel building step in reciprocal space.

As shown in Table III, the wall time of the computing MGP kernel is about 200 times larger than a single call for the evaluation of the potential. Fortunately, the kernel is evaluated only once at the beginning of the simulation. We recall that the computational cost in the initial MGP kernel is strictly linear scaling with system size instead of quasilinear [i.e., N ln(N) scaling with N being the number of real-space grid points] in the evaluation of the potential. As a result, for small systems, the total wall time for optimizing electron density with MGP KEDFs is always between WT and WGC’s wall times. Increasing the number of grid points (or the size of the system), the wall time should approach WT (although this should be system dependent).

Among WGC, MGP, and WT KEDFs, WGC is the most expensive KEDF in terms of total wall time. Additionally, we find that MGP KEDFs converge in all cases considered with 12 iterations of a truncated Newton (TN) minimization. By contrast, WGC KEDF is numerically unstable in some cases (this is a well-known limitation of this functional^10,22).

In summary, the computational efficiency of MGP closely resembles WT’s, and both MGP and WT are more computationally efficient than WGC.

We have formulated and implemented MGP, a new nonlocal kinetic energy functional with a density independent kernel. In MGP, the inverse response function of the FEG is functional-integrated to yield a new kernel. In a second step, the kernel is augmented by a “Kinetic electron” which is opposite to the exchange hole.

Our pilot calculations show that MGP improves dramatically over currently available nonlocal functionals with density independent kernels and performs even better than the best available functionals (featuring density dependent kernels) for both metallic and semiconducting phases of Si as well as the ZB phases of nine common III-V semiconductors. Although the results presented here are quite encouraging, we should take them with a grain of salt. In this work, we have only considered bulk systems, and the lingering question is, can MGP deliver similar quality results for finite systems or systems with vacancies? Initial tests of MGP applied to finite systems (isolated clusters) are quite encouraging and will be the subject of a follow-up study.

In light of recent debates criticizing the use of the total electronic energy as the only descriptor relevant for optimizing a density functional,⁵⁴ we have also inspected the ability of MGP to reproduce electron density distributions. We find that MGP performs quite remarkably in this respect. For example, MGP’s densities are much closer to KS-DFT in bonding regions for CD silicon and ZB GaAs compared to the current state-of-the-art HC functional. Interestingly, we find that although for semiconductors and nonsimple metals, WT is unable to reproduce observables related to the energy; it delivers quite accurate electron densities. As MGP is derived from WT, this is probably the reason why MGP delivers accurate densities.

Finally, the benchmarks for computational cost indicate that our new KEDF is almost as computationally efficient as WT which is of similar scaling as GGA functionals.

Our analysis shows that there is room for improvement, particularly for two aspects of the MGP functional. First, the Kinetic electron can be better parametrized following prescriptions that have been very successful in formulating exchange functionals in real space.^38,55 Second, there exist formulations of “jellium with gap” models for the Lindhard function⁵⁶ which have been recently (and successfully) applied to the formulation of GGA KEDF.⁷ These can be extended to the line-integral formulation of MGP. Third, a density-dependent kernel version of MGP can be constructed to satisfy additional exact conditions, such as density and coordinate scaling relations.³² 

See supplementary material for an analysis of different ways of symmetrizing the kernel with respect to the spatial coordinates and the respective electron density dependence.

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1553993. We thank Emily Carter and her group members for insightful discussions.