Next Article in Journal
Hierarchical Epidemic Model on Structured Population: Diffusion Patterns and Control Policies
Next Article in Special Issue
Periodic DFTB for Supported Clusters: Implementation and Application on Benzene Dimers Deposited on Graphene
Previous Article in Journal
Science Mapping the Academic Knowledge on Business Improvement Districts
Previous Article in Special Issue
Accuracy and Precision in Electronic Structure Computation: Wien2k and FPLO
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Kinetic Energy Density Functionals Based on a Generalized Screened Coulomb Potential: Linear Response and Future Perspectives

by
Eduardo Fabiano
1,2,
Fulvio Sarcinella
2,3,
Lucian A. Constantin
4 and
Fabio Della Sala
1,2,*
1
Institute for Microelectronics and Microsystems (CNR-IMM), Via Monteroni, Campus Unisalento, 73100 Lecce, Italy
2
Center for Biomolecular Nanotechnologies@UNILE, Istituto Italiano di Tecnologia, Via Barsanti 14, 73010 Arnesano, Italy
3
Department of Mathematics and Physics “E. De Giorgi”, University of Salento, Via Arnesano, 73100 Lecce, Italy
4
Istituto di Nanoscienze, Consiglio Nazionale delle Ricerche (CNR-NANO), 41125 Modena, Italy
*
Author to whom correspondence should be addressed.
Computation 2022, 10(2), 30; https://doi.org/10.3390/computation10020030
Submission received: 22 January 2022 / Accepted: 10 February 2022 / Published: 15 February 2022

Abstract

:
We consider kinetic energy functionals that depend, beside the usual semilocal quantities (density, gradient, Laplacian of the density), on a generalized Yukawa potential, that is the screened Coulomb potential of the density raised to some power. These functionals, named Yukawa generalized gradient approximations (yGGA), are potentially efficient real-space semilocal methods that include significant non-local effects and can describe different important exact properties of the kinetic energy. In this work, we focus in particular on the linear response behavior for the homogeneous electron gas (HEG). We show that such functionals are able to reproduce the exact Lindhard function behavior with a very good accuracy, outperforming all other semilocal kinetic functionals. These theoretical advances allow us to perform a detailed analysis of a special class of yGGAs, namely the linear yGGA functionals. Thus, we show how the present approach can generalize the yGGA functionals improving the HEG linear behavior and leading to an extended formula for the kinetic functional. Moreover, testing on several jellium cluster model systems allows highlighting advantages and limitations of the linear yGGA functionals and future perspectives for the development of yGGA kinetic functionals.

1. Introduction

The non-interacting kinetic energy (KE) functional is one of the main quantities of interest in density functional theory [1,2]. Its exact formal definition is readily obtained, following the Levy constrained search formalism [3], as
T s [ n ] = min Ψ n Φ | 1 2 2 | Φ ,
where n is the electron density and Φ is a Kohn-Sham (KS) [4] Slater determinant yielding the density n. This formula allows us to study different important properties of the KE functional and provides an explicit expression for T s in terms of KS orbitals. However, it does not allow us to obtain an explicit expression in terms of the electron density. Therefore, the quest for the KE density functional is still open, also considering the importance of this quantity in many contexts including orbital-free density functional theory (OF-DFT) [5,6,7,8], subsystem DFT [9,10,11,12,13,14,15], and quantum hydrodynamic theory [16,17,18,19]. In addition, semilocal KE functionals have been used in meta-GGA exchange-correlation functionals to remove their orbital dependence [20,21,22,23,24].
Numerous investigations have been dedicated in the last decades to the study of KE functionals [25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72]. Nevertheless, accurate approximations of T s [ n ] are hard to obtain because this quantity usually gives a dominant contribution to the ground-state energy [3] and because of the highly non-local nature of the KE functional [5,51,73,74,75,76,77,78]. For this reason, more recently, machine-learning methods have also been used to develop KE functionals [79,80,81,82,83,84,85,86,87].
A KE functional approximation can be written
T s [ n ] = τ [ n ] ( r ) d r ,
where τ ( r ) is the KE density and, actually, there are two main strategies to approximate τ ( r ) . The simplest one considers the KE density to be a semilocal functional of the density, that is
τ [ n ] ( r ) = τ s e m i l o c a l ( n ( r ) , n ( r ) , 2 n ( r ) , ) .
This approach, which traces back to the nearsightedness principle [88], is computationally efficient because of the local nature of τ ( r ) . However, since it is not explicitly including non-local effects, that are quite relevant for KE functionals, semilocal functionals face several limitations [46,89].
To overcome this problem, the other popular approach used to describe the KE density makes explicit use of a non-local ansatz [47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63]
τ n o n l o c a l ( r ) = n α ( r ) K ( r , r ) n β ( r ) d r ,
where K ( r , r ) is a proper non-local kernel. The presence of the non-local kernel strongly increases the computational cost of the method and raises several practical difficulties. On the other hand, non-local KE functionals are much more accurate [47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63].
In particular, the kernel can be designed in order to reproduce the correct linear response behavior of the homogeneous electron gas (HEG), which has been shown to be a very important property for the KE functional [47,60]. In fact, the KE functional and the linear response function χ have a close relation, given by the equation [5]
F δ 2 T s [ n ] δ n ( r ) δ n ( r ) = 1 χ ( k ) ,
where F denotes the Fourier transform and k is the momentum vector. For the HEG, the linear response function can be computed analytically [5] and is related to the Lindhard function [90]. Hence,
χ H E G ( k ) = ( 3 π 2 n ) 1 / 3 π 2 F L i n d 1 ( η ) ,
F L i n d 1 ( η ) = 1 2 + 1 η 2 4 η ln 1 + η 1 η ,
with η = | k | / [ 2 ( 3 π 2 n ) 1 / 3 ] being a dimensionless momentum. The Lindhard function cannot be accurately mimicked by any semilocal functional, because these all have a polynomial Fourier transform [71]. On the other hand, it is the main property for most of the non-local KE functionals.
Recently, a new class of KE functionals [71] has been proposed to join the advantages of the semilocal methods and the good features of the non-local functionals. These functionals, named Yukawa generalized gradient approximation (yGGA), use as input ingredients, beside the density and its gradients, a Yukawa potential
u α ( r ) = n ( r ) e α k F ( r ) | r r | | r r | d r ,
i.e., a screened Coulomb potential with α k F ( r ) = α ( 3 π 2 ) 1 / 3 n ( r ) 1 / 3 as the screening length. In this way, it is possible to include efficiently non-local effects and improve the description of the HEG response function without resorting to the reciprocal space [71] (as it is instead accustomed in non-local functionals).
In this work, we will take a step further in this direction and we will consider a modification of the basic input ingredient of the functional, allowing the Yukawa potential to be computed for a power of the density, i.e.,
u α β ( r ) = n β ( r ) e α k F ( r ) | r r | | r r | d r .
The computational cost of Equation (9) is, at least for the spherical systems considered in this work, the same as the one of the conventional Yukawa potential in Equation (8); thus, it is interesting to investigate if and how the β parameter will impact the accuracy of the linear response and of the resulting KE functional.
Therefore, we will consider yGGA functionals of the general form
τ [ n ] ( r ) = τ ( n ( r ) , n ( r ) , 2 n ( r ) , u α β ( r ) ) .
For these functionals, we will provide a full analytical derivation of the linear response function and we will consider the exact constraints required to reproduce the Lindhard behavior. Finally, we will analyze the role of the β parameter for the description of jellium spheres and we will give insights for the further development of yGGA functionals.

2. Theory

We consider a KE density of the form
τ = C F n 5 / 3 F s ( p , q , y α β ) ,
where F s is the KE enhancement factor and
C F = 3 10 k 0 2 with k 0 = ( 3 π 2 ) 1 / 3 ,
p ( r ) = | n ( r ) | 2 4 k 0 2 n ( r ) 8 / 3 ,
q ( r ) = 2 n ( r ) 4 k 0 2 n ( r ) 5 / 3 ,
y α β ( r ) = 3 π α 2 4 k 0 n ( r ) β 2 / 3 n β ( r ) e α k 0 n ( r ) 1 / 3 | r r | | r r | d r .
The quantity y α β is the main ingredient of the yGGA functional, whereas p (the reduced gradient) and q (the reduced Laplacian) are the conventional ingredients of meta-GGA functionals. The ingredient y α β is proportional to the potential u α β , where the normalization constant has been choosen so that y α β is adimensional and invariant under uniform scaling, as shown in Ref. [71]. Moreover, in the case of a large number of electrons, where the Thomas–Fermi (TF) limit is exact, we have also y α β 1 : thus, the enhancement factor F s must be 1 when p = 0 ,   q = 0 , and y α β = 1 .
While the case of β = 1 has been deeply investigated in Ref. [71], when β 1 , the properties of y α β change. In particular, in the tail of finite systems we have:
y α β 3 π α 2 4 k 0 Q β n β 2 / 3 r ,
Q β = n β ( r ) d r ,
(for β = 1 , Q β is just the number of electrons). Thus, y α β diverges for β > 2 / 3 (as the density vanishes exponentially), whereas for β 2 / 3 , it vanishes.

2.1. Kinetic Energy Potential

Following the same derivation as in Ref. [71], the KE potential of a general modified yGGA functional is
δ T s yGGA δ n ( r ) = v k M G G A ( r ) + v k yGGA , 1 ( r ) + v k yGGA , 2 ( r ) ,
with
v k M G G A ( r ) = τ n ( r ) · τ n + 2 τ 2 n ,
v k yGGA , 1 ( r ) = β n ( r ) β 1 τ u α β ( r ) e α k F ( r ) | r r | | r r | d 3 r ,
v k yGGA , 2 ( r ) = α k F ( r ) 3 n ( r ) τ u α β ( r ) n β ( r ) e α k F ( r ) | r r | d 3 r .
The potential can also be expressed as a function of the enhancement factor using:
d τ d u α β = 9 π α 2 40 k F ( r ) n ( r ) 2 β d F s d y α β .
Note that the potential in Equation (20) is a non-local function of the expression in Equation (22). This means, for example, that a divergence of d τ d u α β in the tail of a finite system will have an impact everywhere in the space. Thus, the enhancement factor F s must be properly defined in all points.
In the limit of a large number of electrons ( k F ), we have
v k yGGA , 1 ( r ) = β n ( r ) β 1 9 π α 2 40 k F ( r ) n ( r ) 2 β d F s d y α β 4 π α 2 k F ( r ) 2 = β 3 10 k F ( r ) 2 d F s d y α β ( r ) ,
v k yGGA , 2 ( r ) = α k F ( r ) 3 n ( r ) 9 π α 2 40 k F ( r ) n ( r ) 2 β d F s d y α β n β ( r ) 8 π α 3 k F ( r ) 3 = 2 10 k F ( r ) 2 d F s d y α β ( r ) ,
which are stable expressions for all values of β . In this case, we also have that the MGGA term of Equation (19) reduces to the first LDA term only. Hence, for the total KE potential, we find
T s yGGA n ( r ) = 5 10 k F ( r ) 2 F s ( y α β ( r ) ) ( 3 10 β 2 10 ) k F ( r ) 2 y α β ( r ) d F s d y α β ( r ) .
With the condition F s 1 , we have that the total KE potential of Equation (25) reduces to τ / n = ( 5 / 10 ) k F ( r ) = v T F ( r ) , where v T F is the Thomas–Fermi potential. Thus, the condition that in the TF limit F s y α β 1 yields also a correct potential.

2.2. Linear Response of a yGGA Functional

Following Refs. [71,91], we compute the linear response considering the perturbed density n = n 0 + n k e i k · r . Note that the derivation in Ref. [71] was limited to a specific class of yGGA functionals (the linear yGGA functional, see Section 3 in the following) and with β = 1 . Hereafter, a completely general derivation is presented.
Expanding the resulting perturbed KE density in power series, we obtain the linear response χ = 1 / F as twice the coefficient of the inverse of the second-order term. For simplicity, we can evaluate the whole expression in r = 0 , since the HEG is homogeneous and isotropic. Then, we have n = n 0 + n k , | n | 2 = n k 2 k 2 and 2 n = n k k 2 . Therefore, we need to consider
τ = C F ( n 0 + n k ) 5 / 3 F s ( p , q , y α β ) ,
with
p = n k 2 k 2 4 k 0 2 ( n 0 + n k ) 8 / 3 ,
q = n k k 2 4 k 0 2 ( n 0 + n k ) 5 / 3 ,
y α β 1 β k 2 A 2 + k 2 n k n 0 + + β ( β 1 ) 2 A 2 A 2 + 4 k 2 + 2 β A 2 k 2 3 β 2 A 2 ( A 2 + k 2 ) 3 ( A 2 + k 2 ) 2 + β ( β + 1 ) 2 n k n 0 2 ,
where we have set A = α k 0 n 0 1 / 3 . For the derivation of the expression for y α β , see Appendix A.
To proceed, we can expand in powers of n k both factors in Equation (26). Thus, using the notation n k = / n k | n k = 0 , we can write
τ C F n 0 5 / 3 + 5 3 n 0 2 / 3 n k + 5 9 n 0 1 / 3 n k 2 + × × F s ( 0 , 0 , 1 ) + n k F s n k + 1 2 n k 2 F s n k 2 + = = C F n 0 5 / 3 F s ( 0 , 0 , 1 ) + 5 C F n 0 2 / 3 3 F s ( 0 , 0 , 1 ) + C F n 0 5 / 3 n k F s n k + + 5 C F n 0 1 / 3 9 F s ( 0 , 0 , 1 ) + 5 C F n 0 2 / 3 3 n k F s + C F n 0 5 / 3 2 n k 2 F s n k 2 + O ( n k 3 ) .
The linear response χ = 1 / F is given by twice the coefficient of the inverse of the second-order term. Therefore,
F = 2 5 C F n 0 1 / 3 9 F s ( 0 , 0 , 1 ) + 5 C F n 0 2 / 3 3 n k F s + C F n 0 5 / 3 2 n k 2 F s = = k 0 2 n 0 2 / 3 F s ( 0 , 0 , 1 ) 3 n 0 + n k F s + 3 n 0 10 n k 2 F s .
The corresponding Thomas–Fermi-renormalized linear response is F ¯ = k F F / π 2 = k 0 n 1 / 3 F / π 2 . Then,
F ¯ = 3 n 0 F s ( 0 , 0 , 1 ) 3 n 0 + n k F s + 3 n 0 10 n k 2 F s .
For simplicity of notation in the following, we neglect the subscript α β in the y α β ingredient.
To obtain a more explicit expression for Equation (32), we use the chain rule
n k F s = D p n k p + D q n k q + D y n k y ,
where we have employed the notation ( i = p , q , y )
D i F s i | ( p , q , y ) = ( 0 , 0 , 1 ) .
Hence, substituting the values for n k p , n k q , and n k y , we find
n k F s = D q k 2 4 k o 2 n 0 5 / 3 D y 1 n 0 β k 2 A 2 + k 2 .
For the second derivative, we have
n k 2 F s = n k p n k F s p + q n k F s q + y n k F s y = = n k 2 p D p + n k 2 q D q + n k 2 y D y + n k p 2 D p p + n k q 2 D q q + n k y 2 D y y + + 2 n k p n k q D p q + 2 n k p n k y D p y + 2 n k q n k y D q y .
Since n k p = 0 , this immediately simplifies to
n k 2 F s = n k 2 p D p + n k 2 q D q + n k 2 y D y + n k q 2 D q q + n k y 2 D y y + 2 n k q n k y D q y .
Substituting the values for the various derivatives, we find
n k 2 F s = k 2 2 k 0 2 n 0 8 / 3 D p + 5 k 2 6 k 0 2 n 0 8 / 3 D q + + 2 n 0 2 β ( β 1 ) 2 A 2 A 2 + 4 k 2 + 2 β A 2 k 2 3 β 2 A 2 ( A 2 + k 2 ) 3 ( A 2 + k 2 ) 2 + β ( β + 1 ) 2 D y + + k 4 16 k 0 4 n 0 10 / 3 D q q + 1 n 0 2 β 2 k 4 ( A 2 + k 2 ) 2 D y y + β k 4 2 k 0 2 n 0 8 / 3 ( A 2 + k 2 ) D q y .
Finally, using Equations (35) and (38) into Equation (32), we obtain the formula
F ¯ = F s ( 0 , 0 , 1 ) + 9 20 k 2 k 0 2 n 0 2 / 3 D p + 9 160 k 4 n 0 4 / 3 k 0 4 D q q + + 3 5 3 β ( β 1 ) 2 A 2 A 2 + 4 k 2 3 β A 2 ( β A 2 + k 2 ) + b k 2 ( 3 β A 2 + 5 k 2 ) ( A 2 + k 2 ) 2 + 3 β ( β + 1 ) 2 D y + + 9 10 β 2 k 4 ( A 2 + k 2 ) 2 D y y + 9 20 β k 4 k 0 2 n 0 2 / 3 ( A 2 + k 2 ) D q y ,
which, with the substitution η = α k / ( 2 A ) and after some algebra, assumes the final form
F ¯ = F s ( 0 , 0 , 1 ) + 9 5 η 2 D p + 9 10 η 4 D q q + 36 5 β η 4 α 2 + 4 η 2 D q y + 9 10 β ( β 1 ) 1 + α 2 α 2 + 16 η 2 2 α 2 α 2 + 4 η 2 D y + + 24 5 β η 4 α 2 + 4 η 2 2 3 β D y y 4 D y .
We note that this expression is well defined and continuous for any value of α and β . However, the case β = 1 is a special one. In fact, in this case, F ¯ loses its dependence on D y , being dependent only on the linear combination 3 β D y y 4 D y as well as on F s ( 0 , 0 , 1 ) , D p , D q q , D q y (here, we only consider the dependence on the degrees of freedom related to the modeling of the enhancement factor; α and β , which are related to the definition of y, are considered as additional parameters). Therefore, when β = 1 , there is a reduced parametric flexibility to optimize the linear response function through the modeling of the enhancement factor. As we will see, this has consequences for the possibility to impose the correct asymptotic behavior to F ¯ .

Asymptotic Behavior

The asymptotic expansions of the exact response function are [5]
F ¯ L i n d = 1 + η 2 3 + 8 45 η 4 + O ( η 6 ) η 0 ,
F ¯ L i n d = 3 η 2 3 5 + O ( η 2 ) η .
From Equation (40), we get for the small- η limit
F ¯ = F s ( 0 , 0 , 1 ) + 9 5 D p 4 β ( β 1 ) α 2 D y η 2 + + 3 5 3 2 D q q + 12 β α 2 D q y + 16 α 4 β ( 21 β 23 ) D y + 24 α 4 β 2 D y y η 4 ,
and in the large- η limit
F ¯ = 9 10 D q q η 4 + 9 5 D p + β D q y η 2 + + 3 5 5 3 F s ( 0 , 0 , 1 ) + β ( 3 β 7 ) 2 D y + 3 β 2 2 D y y 3 4 α 2 β D q y .
Therefore, we have the following asymptotic conditions
F s ( 0 , 0 , 1 ) = 1 η 0 O ( η 0 ) ,
D p 4 β ( β 1 ) α 2 D y = 5 27 η 0 O ( η 2 ) ,
3 2 D q q + 12 β α 2 D q y +
+ 16 α 4 β ( 21 β 23 ) D y + 24 α 4 β 2 D y y = 8 27 η 0 O ( η 4 ) ,
9 10 D q q = 0 η O ( η 4 ) ,
9 5 D p + β D q y = 3 η O ( η 2 ) ,
5 3 F s ( 0 , 0 , 1 ) + β ( 3 β 7 ) 2 D y +
+ 3 β 2 2 D y y 3 4 α 2 β D q y = 1 η O ( η 0 ) .
As we saw, in the most general case (i.e., for any value of β ), the linear response function has only five degrees of freedom ( F s ( 0 , 0 , 1 ) , D q q , D p , D q y , 3 β D y y 4 D y ) to satisfy these conditions. Thus, initially, we chose to consider Equations (44)–(48), neglecting for the moment the η 0 order in the large- η limit. With this choice, we obtain
F s ( 0 , 0 , 1 ) = 1 ,
D q q = 0 ,
D p = 5 27 + 4 α 2 β ( β 1 ) D y ,
D q y = 40 27 β 4 α 2 ( β 1 ) D y ,
3 β D y y 4 D y = α 2 ( α 2 60 ) 27 β 36 ( β 1 ) D y .
These conditions fix the asymptotic behavior of the response function up to the η 4 order in the short-range limit and to the η 2 order in the long-range one. Using these values into Equation (40), we find
F ¯ = 34560 η 8 + Q 6 η 6 + Q 4 η 4 + 15 α 4 η 2 α 2 + 72 + 45 α 6 45 α 2 + 4 η 2 2 α 2 + 16 η 2 ,
Q 6 = 16 8 α 4 + 255 α 2 5832 β ( β 1 ) D y + 720 ,
Q 4 = 8 α 2 α 4 + 45 α 2 + 810 .
For large η values, Equation (55) behaves as
F ¯ = 3 η 2 + Δ ,
Δ = ( 1 / 90 ) α 4 ( 81 / 10 ) D y β ( β 1 ) ( 4 / 3 ) α 2 + 1 .
For β = 1 , the exact condition Δ = 3 / 5 can only be obtained for specific values of α , as the dependence from D y vanishes. For β 1 , instead, we can solve for D y obtaining
D y = α 4 120 α 2 + 144 729 β ( β 1 ) .
In this general case, we have also
D y y = 13 α 4 3180 α 2 + 5760 + ( 9 α 4 + 2700 α 2 5184 ) β 2187 β 2 ( β 1 ) .
Equations (60) and (61), together with Equations (50)–(53), fix the asymptotic behavior of the response function up to the η 4 order in the small- η limit and down to the η 0 order in the large- η one.
We note that the expressions for D y and D y y diverge at β = 1 : This reflects the fact that, as discussed above, for β = 1 , it is not always possible to fulfill all the asymptotic conditions. Nevertheless, all equations are well defined for any other value of β as well in the limit β 1 . In fact, the divergence in the denominator in this limit, the contribution D y , becomes negligibly small in Equation (40) as well as in Equations (52) and (53); moreover, the linear combination 3 β D y y 4 D y is always well defined, since it does not diverge for any β > 0 .
Substitution of the asymptotic conditions into Equation (40) finally yields the general yGGA response
F ¯ g e n = 34560 η 8 + 432 η 6 ( 45 α 2 16 ) + 8 α 2 η 4 ( α 4 + 45 α 2 + 810 ) + 15 α 4 η 2 ( α 2 + 72 ) + 45 α 6 45 ( α 2 + 4 η 2 ) 2 ( α 2 + 16 η 2 ) .
Note that, remarkably, this formula displays no dependence on β . Thus, β can be used to optimize the functional beyond the linear response regime. The value of α can instead be optimized by fitting to match as close as possible the exact Lindhard function. To this purpose, we can minimize the quantity
σ = w ( η ) 1 F ¯ l i n d ( η ) 1 F ¯ ( η ) d η ,
where F ¯ l i n d is the normalized Lindhard linear response function and
w ( η ) = e μ ( η 1 ) 2 ,
with μ = 2 , is a weighting function. The use of this weighting function allows to focus the fit in the region close to η = 1 instead of the asymptotic ones that are already included by construction. We remark that the results are only weakly dependent on the choice of μ .
Figure 1 reports the values of the errors with respect to the exact Lindhard function for different values of α . We see that the best reproduction of the Lindhard function is achieved for α = 3.31 , with σ = 0.0513 . However, a quite broad range of values allows us to attain a low error; in particular, for α [ 2.26 , 4.66 ] , the response is always better than the one obtained by the yuk2 functional [71].
The linear response functions for various cases are reported in Figure 2. Inspection of the plot shows that they are all quite similar, as already suggested by the considerations above. For α = 3.31 , we obtain at η = 1 , where the Lindhard function shows a derivative singularity, that 1 / F ¯ = 0.47 , which is very close to the exact value 1 / F ¯ = 1 / 2 [5]. Actually, it is also possible also to satisfy 1 / F ¯ = 1 / 2 exactly using α = 3.64 , even if the global accuracy is somehow smaller in this case ( σ = 0.0527 ).

3. Linear yGGA Functionals

We consider the simplest case of yGGA functional by taking functionals that have an enhancement factor with the general form
F s linyGGA ( p , q , y ) = 1 G 0 + 5 3 p + y ( G 0 + G ( p , q ) ) .
i.e., a linear dependence on y. Note that these are closely related with the yGGAs defined in Ref. [71], which are recovered if G 0 = 1 and β = 1 . However, in this work, we lift this restriction.
From Equation (65), we have that F s ( 0 , 0 , 1 ) = 1 is satisfied by construction if G ( 0 , 0 ) = 0 and D q q = G q q , D p = 5 / 3 + G p , D q y = G q , D y y = 0 , D y = G 0 , where we have used the short-hand notation G i = G ( 0 , 0 , 1 ) / i , and G i j = G ( 0 , 0 , 1 ) / i j . Using Equations (50)–(53), we then obtain
G q q = 0
G p = 40 27 + 4 α 2 β ( β 1 ) G 0 = β G q .
The simplest functional satisfying these conditions is the one with
G ( p , q ) = 40 27 β 4 α 2 ( β 1 ) G 0 ( q β p ) x
G 0 = α 2 ( α 2 60 ) 108 β 9 β 10 .
Equation (68) has been derived from Equation (54) using the fact that, for functionals defined by Equation (65), D y y = 0 . The Equations (65), (68), and (69) define the yuk2 β functional, which reproduces the Lindhard functional at small η up to fourth order and for large η behaves as
F ¯ y u k 2 β 3 η 2 + Δ y u k 2 β
Δ y u k 2 β = 1 + ( 2700 β + 3180 ) α 2 360 ( 9 β 10 ) + ( 9 β 13 ) α 4 360 ( 9 β 10 ) .
Although the exact value Δ y u k 2 β = 3 / 5 can be obtained for specific values of α and β , this condition does not have a big impact on the overall accuracy of the linear response. We consider the general error indicator of the linear response of yuk2 β with respect to the Lindhard function (Equation (63)). Results for different values of the α and β are reported as a colormap in Figure 3.
We find that, as seen in the previous section, there is a quite broad range of α values that yield small errors (blue areas in the plot). Moreover, the errors are almost independent on β , except for values β 1.1 , where a discontinuity in the linear response behavior occurs.
Moreover, from the results of the previous section, we can fix α = 3.31 , such that G 0 4.484 / [ β ( 10 9 β ) ] , which is shown in Figure 4.
Note that G 0 is positive up to β = 10 / 9 , where it has a pole, and it has a minimum at β = 5 / 9 . We remark that instead setting G 0 = 1 and β = 1 , then Equation (68) readily yields α = 1.36 . Thus, the yuk2 β functional immediately reduces to the yuk2 functional of Ref. [71].
The yuk2 β functional is just a simple ansatz recovering an accurate linear response behavior. However, because it uses a linear dependence on p and q, it may display severe drawbacks in real applications. In particular, the positivity of the Pauli KE density must be ensured [38,92], which is not the case for Equation (65). In fact, the quantity
1 + w = 1 G 0 + y ( G 0 + x )
is not always positive. Thus, we define the yuk3 β functional as
F s yuk 3 β = 5 3 p + T w ,
with T ( w ) being a positive function such that for w 0 , we have T ( w ) = 1 + w + O ( x 3 ) , such as the one considered in Ref. [71]. Thus, lacking any quadratic term, the yuk3 β has the same linear response of yuk2 β . In this way, although the functional F s yuk 3 β is not truly a linear function of y, it practically behaves as a linear function of y because D y y = 0 . True non-linear yGGA functionals, with D y y 0 , are much more complicated, and they will be considered elsewhere.

4. Computational Details

Densities and enhancement factors were computed in post-SCF fashion using the orbitals of self-consistent Kohn–Sham calculations. All Kohn–Sham calculations have been performed using an in-house developed code solving numerically the spherical symmetry Kohn–Sham problem with the local density approximation being used for the exchange-correlation functional.
Neutral jellium spheres with N = 40, 92, 254 electrons and Wigner–Seitz radius r s = 2, 3, 4, 5, and 6 were considered. They are characterized by the positive background density
n + ( r ) = 3 / ( 4 π r s 3 ) , r < R 0 , r R .
The post-processing of the Kohn–Sham orbitals to obtain the quantities G ( p , q ) , y, and w has been carried out by an additional in-house software which evaluates the required derivatives and integrals in real space using the same grid as the one used in the Kohn–Sham calculations.

5. Results

We considered different pairs of α and β values for the y ingredient, as reported in Table 1 (see also Figure 3). The first pair ( α = 1.36 , β = 1 ) is the one of the original yuk2 functional of Ref. [71]. The second, third, and fourth pairs consider the α value suggested by Figure 1 and several values of β . The last pair considers the minimizing β value of Figure 4 and a corresponding α such that G 0 = 1 . All these pairs give a quite accurate description of the HEG linear response: the highest accuracy is obtained for ( α , β ) = (3.31, 2) and (2.34, 5/9), as reported in the last column of Table 1.
In Figure 5, we report, in panels (a), (b), (c) and (d) respectively, the electronic density, the function x defined in Equation (68), as well as the values of y and w for a jellium sphere with N = 254 electrons and r s = 4 .
Looking at Figure 5, we see in panel (b) that for β 1 , the values of x are negative in the tail, as also discussed in Appendix C. Moreover, panel (c) shows that the ingredient y is close to 1 inside the jellium spheres for all ( α , β ) pairs but in the tail, it diverges for β 1 . On the other hand, for β = 2 / 3 , it reaches a constant (but it will slowly vanish in very far regions), while for β 2 / 3 , it soon decreases to zero, as discussed in Section 2. Finally, panel (d) shows the ingredient w. Inside the sphere, the values of w are in the range | w | < 0.3 for all the ( α , β ) combinations. The behavior in the tail follows the one of x but for β = 1 and α = 3.31 , there is a peak before the negative divergence.
These results indicate that for the total kinetic energy, which is mainly influenced by the functional behavior inside the jellium sphere, all the cases considered here can be expected to yield similar results, which are in line with the performance of the yuk3 functional [71]. On the other hand, for the potential, we have a contribution F s / y F s / w · w / y . However, this latter term for functionals with β < 1 is largely divergent ( w ( r ) is diverging but y ( r ) is not); similarly, it diverges in the ( α = 3.31 , β = 1 ) case, where w strongly oscillates around r = 30 bohr. Thus, we can expect, for these values of the parameters, major problems on the KE potential.
To understand better the role of β on the KE functional, we can perform a further analysis of the “exact” Pauli KE enhancement factor
F p exact τ p K S ( r ) / τ T F ( r ) ,
where τ p K S is the Pauli KE density corresponding to the Kohn–Sham positive-defined KE density and τ T F is the Thomas–Fermi density. Hence, we have computed F p exact for several jellium clusters with different numbers of electrons ( N = 40 , 92 , 254 ) and r s = 2, 3, 4, 5, and 6 on all grid points, and we have plotted it versus the corresponding values of w. In fact, Equation (73) implies that that the Pauli enhancement factor can be written as a universal single-valued function of w. However, this is an ansatz, which must be verified; if it is correct, then we must obtain a unique value of F p for each w, also for different systems.
Figure 6 shows the plot of F p exact vs. w for the various pairs of parameters considered in this work. It is evident that the cases with β < 1 and the case with ( α = 3.31 , β = 1 ) do not yield single-valued functions but rather display multiple branches of F p exact for all values of w. This confirms that using these parameters, it is not possible to obtain an accurate description of the kinetic functional. On the contrary, the yuk2 parameterization ( α = 1.35 , β = 1 ) and the one with ( α = 3.31 , β = 2 ) show a nice single line for w 0 . For negative w values, they display a multi-valued behavior; however, for these cases, the w < 0 region corresponds to the density tail region. Therefore, the effect, at least on the computation of the energy, is minor. Thus, we can infer that both functionals will be accurate for the kinetic energies (this is indeed the case for yuk2 [71]) but may yield some oscillations in the asymptotic part of the potential because of the multi-valued behavior at w < 0 . This latter feature is possibly significantly reduced for ( α = 3.31 , β = 2 ) that displays the lowest grade of multi-valuedness at negative w values.
In a future work, we will develop an approach that will allow us to obtain an analytical F ( w ) function in order to build up an accurate KE functional and KE potential. However, the results reported in Figure 6 also indicate that the ansatz in Equation (73) may be not accurate enough and further studies are required, in particular, to investigate the effect of the non-linear term ( y 2 ), which is already important for the HEG linear-response and can be even more important for the description of finite systems.

6. Conclusions

In this paper, we have investigated non-interacting kinetic energy functionals depending on a generalized Yukawa potential, i.e., a screened Coulomb potential of the electronic density raised to the power of β . The use of this input ingredient allows us to introduce in an efficient way non-local features into the functional.
In particular, we have derived the exact homogeneous electron gas linear response behavior of generalized yGGA functionals and, by comparing to the Lindhard function, we have derived exact asymptotic constraints for the functional.
In particular, it turned out that β = 1 is a very special case, and the Lindhard asymptotic constraints can only be satisfied for a specific value of the screening parameter ( α ). Moreover, the final linear response of yGGA functionals satisfying the low and high-wavevector Lindhard properties does not depend on β , which can be then used as an additional degree of freedom to model systems beyond the linear-response regime.
We have used the developed theory to extend the work reported in Ref. [71] and investigate in detail the simplest class of yGGA functionals, namely the linear yGGAs, i.e., those yGGA functionals depending only linearly on y. We have found that although this class of functionals can satisfy rather accurately the linear response constraints, it is not flexible enough to perform very accurately for both the kinetic energy and potential computation. This is mainly due to the fact that imposing the linear response behavior implies, in these simple functionals, that the Pauli enhancement factor is a function of a well-defined combination of the density ingredients (the w ingredient of Equation (72)), but this is not sufficient to describe the non-local nature of the Pauli kinetic energy. Although for some wise choices of the parameters, this effect can be minimized; this is an intrinsic limitation of the linear yGGA functional class.
Therefore, future work will focus on the development of more sophisticated functionals forms using an explicit non-linear dependence on the Yukawa potential. The theoretical framework established in this paper will possibly allow us to develop more efficient and broadly applicable kinetic energy functionals.

Author Contributions

Conceptualization and methodology, E.F. and F.D.S.; software, F.S.; investigation, E.F., F.S. and F.D.S.; writing, E.F., F.S., L.A.C. and F.D.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Weak Perturbation of y

To compute the perturbation of y defined in Equation (15), we first consider the generalized Yukawa potential
u α ( r ) = n β ( r ) e A ( r ) | r r | | r r | d r ,
where we have set A = α k 0 n 0 1 / 3 . Upon the density perturbation, this becomes
u α = n 0 β 1 + e i k · r n k n 0 β e A 1 + n k n 0 1 / 3 r r d r .
Expanding in series and keeping only terms up to second order, we find
u α n 0 β I 1 + β I 2 1 3 A I 3 n k n 0 + β ( β 1 ) 2 I 4 + A 18 ( A I 5 + 2 I 3 ) A β 3 I 6 n k n 0 2 ,
with
I 1 e A r r d r
I 2 e i k · r A r r d r
I 3 e A r d r
I 4 e 2 i k · r A r r d r
I 5 r e A r d r
I 6 e i k · r A r d r .
Substituting the values of the integrals (see Appendix B), we find
u α 4 π n 0 β A 2 + 4 π n 0 β β A 2 + k 2 2 3 A 2 n k n 0 + + 4 π n 0 β β ( β 1 ) 2 1 A 2 + 4 k 2 + 5 9 A 2 2 A 2 β 3 ( A 2 + k 2 ) 2 n k n 0 2 .
Now, we can consider the perturbation of y, that is
y = 3 π α 2 4 k 0 ( n 0 + n k ) β 2 / 3 u α [ n 0 + n k ] ( r ) = = 3 π α 2 n 0 2 / 3 β 4 k 0 1 + n k n 0 2 / 3 β u α [ n 0 + n k ] ( r ) .
Using the result of Equation (A9) and expanding to second-order, we find
y 1 β k 2 A 2 + k 2 n k n 0 + + β ( β 1 ) 2 A 2 A 2 + 4 k 2 + 2 β A 2 k 2 3 β 2 A 2 ( A 2 + k 2 ) 3 ( A 2 + k 2 ) 2 + β ( β + 1 ) 2 n k n 0 2 .

Appendix B. Integrals for the Generalized Yukawa Expansion

For integral I 1 , I 3 and I 5 , we use that
E [ A , n ] = 0 ( r ) 2 r n e A r d r = Γ ( 3 + n ) A 3 + n .
Thus, I 1 = 4 π E [ A , 1 ] = 4 π A 2 , I 3 = 4 π E [ A , 0 ] = 8 π A 3 , I 5 = 4 π E [ A , 1 ] = 24 π A 4 .
To compute the I 2 , I 4 , and I 6 integrals, we can choose the axis such that k is aligned with the z axis; then, we have that k · r = k r cos θ . Hence, we can write
F [ k , n ] = 0 2 π d ϕ 0 π sin θ d θ 0 d r r 2 ( r ) n e A r e i k r cos θ = = 2 π 0 d r r 2 ( r ) 2 e A r 0 π d θ sin θ e i k r cos θ = = 2 π i k 0 ( r ) n r e A r e i k r e i k r d r = 2 π i k ( E [ A i k , n ] E [ A + i k , n ] ) = 4 π k [ E [ A i k , n 1 ] ] .
Thus, I 2 = F [ k , 1 ] = 4 π A 2 + k 2 and I 6 = F [ k , 0 ] = 8 π A ( A 2 + k 2 ) 2 .
The I 4 integral is similar to the integral I 2 with the substitution k 2 k . Then,
I 4 = 4 π A 2 + 4 k 2 .

Appendix C. Asymptotics

For an exponential spherical density ρ ( r ) = A exp ( 2 Z r ) , where Z = 2 ϵ H , we have
p Z 2 1 k 0 2 1 ρ 2 / 3
q Z ( Z r 1 ) r 1 k 0 2 1 ρ 2 / 3
so that
q p Z r 1 k 0 2 1 ρ 2 / 3 = p r Z .
Thus, q p is negative in the tail; i.e., q is very large in the tail but smaller than p. If we consider q β p with β > 1 , then q β p will be more negative. If we consider q β p with β < 1 , then q β p will be positive.

References

  1. Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864–B871. [Google Scholar] [CrossRef] [Green Version]
  2. Dreizler, R.M.; Gross, E.K.U. Density Functional Theory; Springer: Berlin/Heidelberg, Germany, 1990. [Google Scholar]
  3. Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. Proc. Nat. Acad. Sci. USA 1979, 76, 6062–6065. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Kohn, W.; Sham, L.J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. [Google Scholar] [CrossRef] [Green Version]
  5. Wang, Y.A.; Carter, E.A. Orbital-free kinetic-energy density functional theory. In Theoretical Methods in Condensed Phase Chemistry; Springer: Berlin/Heidelberg, Germany, 2002; pp. 117–184. [Google Scholar]
  6. Wesolowski, T.A.; Wang, Y.A. (Eds.) Recent Progress in Orbital-Free Density Functional Theory; World Scientific: Singapore, 2013. [Google Scholar] [CrossRef] [Green Version]
  7. Gavini, V.; Bhattacharya, K.; Ortiz, M. Quasi-continuum orbital-free density-functional theory: A route to multi-million atom non-periodic DFT calculation. J. Mech. Phys. Sol. 2007, 55, 697–718. [Google Scholar] [CrossRef]
  8. Witt, W.C.; Beatriz, G.; Dieterich, J.M.; Carter, E.A. Orbital-free density functional theory for materials research. J. Mater. Res. 2018, 33, 777–795. [Google Scholar] [CrossRef]
  9. Götz, A.W.; Beyhan, S.M.; Visscher, L. Performance of Kinetic Energy Functionals for Interaction Energies in a Subsystem Formulation of Density Functional Theory. J. Chem. Theory Comput. 2009, 5, 3161–3174. [Google Scholar] [CrossRef]
  10. Jacob, C.R.; Neugebauer, J. Subsystem density-functional theory. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 325–362. [Google Scholar] [CrossRef]
  11. Wesolowski, T.A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891–5928. [Google Scholar] [CrossRef] [Green Version]
  12. Neugebauer, J. Chromophore-specific theoretical spectroscopy: From subsystem density functional theory to mode-specific vibrational spectroscopy. Phys. Rep. 2010, 489, 1–87. [Google Scholar] [CrossRef]
  13. Krishtal, A.; Sinha, D.; Genova, A.; Pavanello, M. Subsystem density-functional theory as an effective tool for modeling ground and excited states, their dynamics and many-body interactions. J. Phys. Condens. Matter 2015, 27, 183202. [Google Scholar] [CrossRef]
  14. Laricchia, S.; Fabiano, E.; Della Sala, F. Frozen density embedding with hybrid functionals. J. Chem. Phys. 2010, 133, 164111. [Google Scholar] [CrossRef] [PubMed]
  15. Śmiga, S.; Fabiano, E.; Laricchia, S.; Constantin, L.A.; Della Sala, F. Subsystem density functional theory with meta-generalized gradient approximation exchange-correlation functionals. J. Chem. Phys. 2015, 142, 154121. [Google Scholar] [CrossRef] [Green Version]
  16. Toscano, G.; Straubel, J.; Kwiatkowski, A.; Rockstuhl, C.; Evers, F.; Xu, H.; Asger Mortensen, N.; Wubs, M. Resonance Shifts and Spill-out Effects in Self-Consistent Hydrodynamic Nanoplasmonics. Nat. Commun. 2015, 6, 7132. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Ciracì, C.; Della Sala, F. Quantum hydrodynamic theory for plasmonics: Impact of the electron density tail. Phys. Rev. B 2016, 93, 205405. [Google Scholar] [CrossRef] [Green Version]
  18. Moldabekov, Z.A.; Bonitz, M.; Ramazanov, T.S. Theoretical Foundations of Quantum Hydrodynamics for Plasmas. Phys. Plasmas 2018, 25, 031903. [Google Scholar] [CrossRef]
  19. Baghramyan, H.M.; Della Sala, F.; Ciracì, C. Laplacian-Level Quantum Hydrodynamic Theory for Plasmonics. Phys. Rev. X 2021, 11, 011049. [Google Scholar] [CrossRef]
  20. Mejia-Rodriguez, D.; Trickey, S.B. Deorbitalization strategies for meta-generalized-gradient-approximation exchange-correlation functionals. Phys. Rev. A 2017, 96, 052512. [Google Scholar] [CrossRef] [Green Version]
  21. Mejia-Rodriguez, D.; Trickey, S.B. Deorbitalized meta-GGA exchange-correlation functionals in solids. Phys. Rev. B 2018, 98, 115161. [Google Scholar] [CrossRef] [Green Version]
  22. Mejía-Rodríguez, D.; Trickey, S.B. Meta-GGA performance in solids at almost GGA cost. Phys. Rev. B 2020, 102, 121109. [Google Scholar] [CrossRef]
  23. Tran, F.; Kovács, P.; Kalantari, L.; Madsen, G.K.H.; Blaha, P. Orbital-free approximations to the kinetic-energy density in exchange-correlation MGGA functionals: Tests on solids. J. Chem. Phys. 2018, 149, 144105. [Google Scholar] [CrossRef] [Green Version]
  24. Tran, F.; Baudesson, G.; Carrete, J.; Madsen, G.K.H.; Blaha, P.; Schwarz, K.; Singh, D.J. Shortcomings of meta-GGA functionals when describing magnetism. Phys. Rev. B 2020, 102, 024407. [Google Scholar] [CrossRef]
  25. Acharya, P.K.; Bartolotti, L.J.; Sears, S.B.; Parr, R.G. An atomic kinetic energy functional with full Weizsacker correction. Proc. Nat. Acad. Sci. USA 1980, 77, 6978–6982. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Thakkar, A.J. Comparison of kinetic-energy density functionals. Phys. Rev. A 1992, 46, 6920. [Google Scholar] [CrossRef]
  27. Lembarki, A.; Chermette, H. Obtaining a gradient–corrected kinetic–energy functional from the Perdew–Wang exchange functional. Phys. Rev. A 1994, 50, 5328. [Google Scholar] [CrossRef]
  28. Tran, F.; Wesołowski, T.A. Link between the kinetic-and exchange-energy functionals in the generalized gradient approximation. Int. J. Quant. Chem. 2002, 89, 441–446. [Google Scholar] [CrossRef]
  29. Constantin, L.A.; Fabiano, E.; Laricchia, S.; Della Sala, F. Semiclassical neutral atom as a reference system in density functional theory. Phys. Rev. Lett. 2011, 106, 186406. [Google Scholar] [CrossRef]
  30. Laricchia, S.; Fabiano, E.; Constantin, L.A.; Della Sala, F. Generalized Gradient Approximations of the Noninteracting Kinetic Energy from the Semiclassical Atom Theory: Rationalization of the Accuracy of the Frozen Density Embedding Theory for Nonbonded Interactions. J. Chem. Theory Comput. 2011, 7, 2439–2451. [Google Scholar] [CrossRef]
  31. Karasiev, V.V.; Chakraborty, D.; Shukruto, O.A.; Trickey, S.B. Nonempirical generalized gradient approximation free-energy functional for orbital-free simulations. Phys. Rev. B 2013, 88, 161108(R). [Google Scholar] [CrossRef] [Green Version]
  32. Borgoo, A.; Tozer, D.J. Density scaling of noninteracting kinetic energy functionals. J. Chem. Theory Comput. 2013, 9, 2250–2255. [Google Scholar] [CrossRef]
  33. Trickey, S.B.; Karasiev, V.V.; Chakraborty, D. Comment on “Single-point kinetic energy density functionals: A pointwise kinetic energy density analysis and numerical convergence investigation”. Phys. Rev. B 2015, 92, 117101. [Google Scholar] [CrossRef] [Green Version]
  34. Xia, J.; Carter, E.A. Reply to “Comment on ‘Single-point kinetic energy density functionals: A pointwise kinetic energy density analysis and numerical convergence investigation’”. Phys. Rev. B 2015, 92, 117102. [Google Scholar] [CrossRef] [Green Version]
  35. Constantin, L.A.; Fabiano, E.; Śmiga, S.; Della Sala, F. Jellium-with-gap model applied to semilocal kinetic functionals. Phys. Rev. B 2017, 95, 115153. [Google Scholar] [CrossRef] [Green Version]
  36. Luo, K.; Karasiev, V.V.; Trickey, S. A simple generalized gradient approximation for the noninteracting kinetic energy density functional. Phys. Rev. B 2018, 98, 041111(R). [Google Scholar] [CrossRef] [Green Version]
  37. Lehtomäki, J.; Lopez-Acevedo, O. Semilocal kinetic energy functionals with parameters from neutral atoms. Phys. Rev. B 2019, 100, 165111. [Google Scholar] [CrossRef] [Green Version]
  38. Perdew, J.P.; Constantin, L.A. Laplacian-level density functionals for the kinetic energy density and exchange-correlation energy. Phys. Rev. B 2007, 75, 155109. [Google Scholar] [CrossRef] [Green Version]
  39. Karasiev, V.V.; Jones, R.S.; Trickey, S.B.; Harris, F.E. Properties of constraint-based single-point approximate kinetic energy functionals. Phys. Rev. B 2009, 80, 245120. [Google Scholar] [CrossRef] [Green Version]
  40. Laricchia, S.; Constantin, L.A.; Fabiano, E.; Della Sala, F. Laplacian-Level kinetic energy approximations based on the fourth-order gradient expansion: Global assessment and application to the subsystem formulation of density functional theory. J. Chem. Theory Comput. 2013, 10, 164–179. [Google Scholar] [CrossRef] [Green Version]
  41. Cancio, A.C.; Stewart, D.; Kuna, A. Visualization and analysis of the Kohn-Sham kinetic energy density and its orbital-free description in molecules. J. Chem. Phys. 2016, 144, 084107. [Google Scholar] [CrossRef] [Green Version]
  42. Cancio, A.C.; Redd, J.J. Visualisation and orbital-free parametrisation of the large-Z scaling of the kinetic energy density of atoms. Mol. Phys. 2017, 115, 618–635. [Google Scholar] [CrossRef] [Green Version]
  43. Constantin, L.A.; Fabiano, E.; Della Sala, F. Semilocal Pauli–Gaussian Kinetic Functionals for Orbital-Free Density Functional Theory Calculations of Solids. J. Phys. Chem. Lett. 2018, 9, 4385–4390. [Google Scholar] [CrossRef] [Green Version]
  44. Constantin, L.A.; Fabiano, E.; Della Sala, F. Performance of Semilocal Kinetic Energy Functionals for Orbital-Free Density Functional Theory. J. Chem. Theory Comput. 2019, 15, 3044–3055. [Google Scholar] [CrossRef] [PubMed]
  45. Śmiga, S.; Constantin, L.A.; Della Sala, F.; Fabiano, E. The Role of the Reduced Laplacian Renormalization in the Kinetic Energy Functional Development. Computation 2019, 7, 65. [Google Scholar] [CrossRef] [Green Version]
  46. Xia, J.; Carter, E.A. Single-point kinetic energy density functionals: A pointwise kinetic energy density analysis and numerical convergence investigation. Phys. Rev. B 2015, 91, 045124. [Google Scholar] [CrossRef] [Green Version]
  47. Wang, L.W.; Teter, M.P. Kinetic-energy functional of the electron density. Phys. Rev. B 1992, 45, 13196. [Google Scholar] [CrossRef]
  48. Perrot, F. Hydrogen-hydrogen interaction in an electron gas. J. Phys. Condens. Matter 1994, 6, 431. [Google Scholar] [CrossRef]
  49. Smargiassi, E.; Madden, P.A. Orbital-free kinetic-energy functionals for first-principles molecular dynamics. Phys. Rev. B 1994, 49, 5220. [Google Scholar] [CrossRef]
  50. García-González, P.; Alvarellos, J.E.; Chacón, E. Nonlocal kinetic-energy-density functionals. Phys. Rev. B 1996, 53, 9509. [Google Scholar] [CrossRef]
  51. García-González, P.; Alvarellos, J.E.; Chacón, E. Kinetic-energy density functional: Atoms and shell structure. Phys. Rev. A 1996, 54, 1897. [Google Scholar] [CrossRef]
  52. García-González, P.; Alvarellos, J.E.; Chacón, E. Nonlocal symmetrized kinetic-energy density functional: Application to simple surfaces. Phys. Rev. B 1998, 57, 4857. [Google Scholar] [CrossRef]
  53. Wang, Y.A.; Govind, N.; Carter, E.A. Orbital-free kinetic-energy functionals for the nearly free electron gas. Phys. Rev. B 1998, 58, 13465. [Google Scholar] [CrossRef] [Green Version]
  54. Wang, Y.A.; Govind, N.; Carter, E.A. Orbital-free kinetic-energy density functionals with a density-dependent kernel. Phys. Rev. B 1999, 60, 16350. [Google Scholar] [CrossRef] [Green Version]
  55. Zhou, B.; Ligneres, V.L.; Carter, E.A. Improving the orbital-free density functional theory description of covalent materials. J. Chem. Phys. 2005, 122, 044103. [Google Scholar] [CrossRef]
  56. Garcia-Aldea, D.; Alvarellos, J. Fully nonlocal kinetic energy density functionals: A proposal and a general assessment for atomic systems. J. Chem. Phys. 2008, 129, 074103. [Google Scholar] [CrossRef]
  57. Garcia-Aldea, D.; Alvarellos, J.E. Approach to kinetic energy density functionals: Nonlocal terms with the structure of the von Weizsäcker functional. Phys. Rev. A 2008, 77, 022502. [Google Scholar] [CrossRef]
  58. Huang, C.; Carter, E.A. Nonlocal orbital-free kinetic energy density functional for semiconductors. Phys. Rev. B 2010, 81, 045206. [Google Scholar] [CrossRef] [Green Version]
  59. Shin, I.; Carter, E.A. Enhanced von Weizsäcker Wang-Govind-Carter kinetic energy density functional for semiconductors. J. Chem. Phys. 2014, 140, 18A531. [Google Scholar] [CrossRef]
  60. Constantin, L.A.; Fabiano, E.; Della Sala, F. Nonlocal kinetic energy functional from the jellium-with-gap model: Applications to orbital-free density functional theory. Phys. Rev. B 2018, 97, 205137. [Google Scholar] [CrossRef] [Green Version]
  61. Mi, W.; Genova, A.; Pavanello, M. Nonlocal kinetic energy functionals by functional integration. J. Chem. Phys. 2018, 148, 184107. [Google Scholar] [CrossRef] [Green Version]
  62. Mi, W.; Pavanello, M. Orbital-free density functional theory correctly models quantum dots when asymptotics, nonlocality, and nonhomogeneity are accounted for. Phys. Rev. B 2019, 100, 041105. [Google Scholar] [CrossRef] [Green Version]
  63. Xu, Q.; Lv, J.; Wang, Y.; Ma, Y. Nonlocal kinetic energy density functionals for isolated systems obtained via local density approximation kernels. Phys. Rev. B 2020, 101, 045110. [Google Scholar] [CrossRef] [Green Version]
  64. Finzel, K. Shell-structure-based functionals for the kinetic energy. Theor. Chem. Acc. 2015, 134, 106. [Google Scholar] [CrossRef]
  65. Finzel, K. Local conditions for the Pauli potential in order to yield self-consistent electron densities exhibiting proper atomic shell structure. J. Chem. Phys. 2016, 144, 034108. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  66. Ludeña, E.V.; Salazar, E.X.; Cornejo, M.H.; Arroyo, D.E.; Karasiev, V.V. The Liu-Parr power series expansion of the Pauli kinetic energy functional with the incorporation of shell-inducing traits: Atoms. Int. J. Quantum Chem. 2018, 118, e25601. [Google Scholar] [CrossRef]
  67. Salazar, E.X.; Guarderas, P.F.; Ludena, E.V.; Cornejo, M.H.; Karasiev, V.V. Study of some simple approximations to the non-interacting kinetic energy functional. Int. J. Quantum Chem. 2016, 116, 1313–1321. [Google Scholar] [CrossRef] [Green Version]
  68. Yao, K.; Parkhill, J. Kinetic Energy of Hydrocarbons as a Function of Electron Density and Convolutional Neural Networks. J. Chem. Theory Comput. 2016, 12, 1139–1147. [Google Scholar] [CrossRef]
  69. Alonso, J.A.; Girifalco, L.A. Nonlocal approximation to the exchange potential and kinetic energy of an inhomogeneous electron gas. Phys. Rev. B 1978, 17, 3735. [Google Scholar] [CrossRef]
  70. Chacón, E.; Alvarellos, J.E.; Tarazona, P. Nonlocal kinetic energy functional for nonhomogeneous electron systems. Phys. Rev. B 1985, 32, 7868. [Google Scholar] [CrossRef]
  71. Sarcinella, F.; Fabiano, E.; Constantin, L.A.; Della Sala, F. Nonlocal kinetic energy functionals in real space using a Yukawa-potential kernel: Properties, linear response, and model functionals. Phys. Rev. B 2021, 103, 155127. [Google Scholar] [CrossRef]
  72. Kumar, S.; Borda, E.L.; Sadigh, B.; Zhu, S.; Hamel, S.; Gallagher, B.; Bulatov, V.; Klepeis, J.; Samanta, A. Accurate parameterization of the kinetic energy functional. J. Chem. Phys. 2022, 156, 024110. [Google Scholar] [CrossRef]
  73. Karasiev, V.V.; Chakraborty, D.; Trickey, S.B. Progress on New Approaches to Old Ideas: Orbital-Free Density Functionals. In Many-Electron Approaches in Physics, Chemistry and Mathematics; Bach, V., Delle Site, L., Eds.; Springer: Cham, Switzerland, 2014; pp. 113–134. [Google Scholar] [CrossRef]
  74. March, N.; Santamaria, R. Non-local relation between kinetic and exchange energy densities in Hartree–Fock theory. Int. J. Quantum Chem. 1991, 39, 585–592. [Google Scholar] [CrossRef]
  75. Della Sala, F.; Fabiano, E.; Constantin, L.A. Kohn-Sham kinetic energy density in the nuclear and asymptotic regions: Deviations from the von Weizsäcker behavior and applications to density functionals. Phys. Rev. B 2015, 91, 035126. [Google Scholar] [CrossRef] [Green Version]
  76. Howard, I.A.; March, N.H.; Van Doren, V.E. r- and p-space electron densities and related kinetic and exchange energies in terms of s states alone for the leading term in the 1/Z expansion for nonrelativistic closed-shell atomic ions. Phys. Rev. A 2001, 63, 062501. [Google Scholar] [CrossRef]
  77. Constantin, L.A.; Fabiano, E.; Della Sala, F. Kinetic and Exchange Energy Densities near the Nucleus. Computation 2016, 4, 19. [Google Scholar] [CrossRef] [Green Version]
  78. Śmiga, S.; Siecińska, S.; Fabiano, E. Methods to generate reference total and Pauli kinetic potentials. Phys. Rev. B 2020, 101, 165144. [Google Scholar] [CrossRef]
  79. Snyder, J.C.; Rupp, M.; Hansen, K.; Müller, K.R.; Burke, K. Finding density functionals with machine learning. Phys. Rev. Lett. 2012, 108, 253002. [Google Scholar] [CrossRef]
  80. Li, L.; Snyder, J.C.; Pelaschier, I.M.; Huang, J.; Niranjan, U.N.; Duncan, P.; Rupp, M.; Müller, K.R.; Burke, K. Understanding machine-learned density functionals. Int. J. Quant. Chem. 2016, 116, 819–833. [Google Scholar] [CrossRef] [Green Version]
  81. Alharbi, F.H.; Kais, S. Kinetic energy density for orbital-free density functional calculations by axiomatic approach. Int. J. Quantum Chem. 2017, 117, e25373. [Google Scholar] [CrossRef]
  82. Seino, J.; Kageyama, R.; Fujinami, M.; Ikabata, Y.; Nakai, H. Semi-local machine-learned kinetic energy density functional with third-order gradients of electron density. J. Chem. Phys. 2018, 148, 241705. [Google Scholar] [CrossRef]
  83. Golub, P.; Manzhos, S. Kinetic energy densities based on the fourth order gradient expansion: Performance in different classes of materials and improvement via machine learning. Phys. Chem. Chem. Phys. 2019, 21, 378–395. [Google Scholar] [CrossRef] [Green Version]
  84. Manzhos, S.; Golub, P. Data-driven kinetic energy density fitting for orbital-free DFT: Linear vs Gaussian process regression. J. Chem. Phys. 2020, 153, 074104. [Google Scholar] [CrossRef]
  85. Meyer, R.; Weichselbaum, M.; Hauser, A.W. Machine Learning Approaches toward Orbital-free Density Functional Theory: Simultaneous Training on the Kinetic Energy Density Functional and Its Functional Derivative. J. Chem. Theory Comput. 2020, 16, 5685–5694. [Google Scholar] [CrossRef] [PubMed]
  86. Imoto, F.; Imada, M.; Oshiyama, A. Order-N orbital-free density-functional calculations with machine learning of functional derivatives for semiconductors and metals. Phys. Rev. Res. 2021, 3, 033198. [Google Scholar] [CrossRef]
  87. Ryczko, K.; Wetzel, S.J.; Melko, R.G.; Tamblyn, I. Toward Orbital-Free Density Functional Theory with Small Data Sets and Deep Learning. J. Chem. Theory Comput. 2022, 18, 1122–1128. [Google Scholar] [CrossRef] [PubMed]
  88. Prodan, E.; Kohn, W. Nearsightedness of electronic matter. Proc. Nat. Acad. Sci. USA 2005, 102, 11635–11638. [Google Scholar] [CrossRef] [Green Version]
  89. Constantin, L.A.; Fabiano, E.; Della Sala, F. Modified Fourth-Order Kinetic Energy Gradient Expansion with Hartree Potential-Dependent Coefficients. J. Chem. Theory Comput. 2017, 13, 4228–4239. [Google Scholar] [CrossRef]
  90. Lindhard, J. On the properties of a gas of charged particles. Dan. Vid. Selsk Mat.-Fys. Medd. 1954, 28, 8. [Google Scholar]
  91. Tao, J.; Perdew, J.P.; Almeida, L.M.; Fiolhais, C.; Kümmel, S. Nonempirical density functionals investigated for jellium: Spin-polarized surfaces, spherical clusters, and bulk linear response. Phys. Rev. B 2008, 77, 245107. [Google Scholar] [CrossRef]
  92. Levy, M.; Ou-Yang, H. Exact properties of the Pauli potential for the square root of the electron density and the kinetic energy functional. Phys. Rev. A 1988, 38, 625–629. [Google Scholar] [CrossRef]
Figure 1. Values of the error function σ (Equation (63)) for the function F ¯ g e n for different values of the α parameter. For comparison, the σ value corresponding to the yuk2 functional, ref. [71] is also reported.
Figure 1. Values of the error function σ (Equation (63)) for the function F ¯ g e n for different values of the α parameter. For comparison, the σ value corresponding to the yuk2 functional, ref. [71] is also reported.
Computation 10 00030 g001
Figure 2. Linear response functions as computed by Equations (55) and (62) [ F ¯ g e n ] for α = 3.31 . For Equation (55), two different choices of β have been considered; for each one, the value of D y has been chosen such to minimize the error σ in Equation (63). The exact Lindhard and the yuk2 responses are also reported for comparison. The inset shows instead the difference between the Lindhard function and the various response functions reported in the plot.
Figure 2. Linear response functions as computed by Equations (55) and (62) [ F ¯ g e n ] for α = 3.31 . For Equation (55), two different choices of β have been considered; for each one, the value of D y has been chosen such to minimize the error σ in Equation (63). The exact Lindhard and the yuk2 responses are also reported for comparison. The inset shows instead the difference between the Lindhard function and the various response functions reported in the plot.
Computation 10 00030 g002
Figure 3. Error with respect to the Lindhard function (Equation (63)) for the yuk2 β functional at various values of the parameters α and β . The black dots denote the positions of the pairs of parameters listed in Table 1.
Figure 3. Error with respect to the Lindhard function (Equation (63)) for the yuk2 β functional at various values of the parameters α and β . The black dots denote the positions of the pairs of parameters listed in Table 1.
Computation 10 00030 g003
Figure 4. Value of G 0 as a function of β for different values of α .
Figure 4. Value of G 0 as a function of β for different values of α .
Computation 10 00030 g004
Figure 5. Electronic density (a) and KE ingredients (bd) for a jellium sphere with N = 254 electrons and r s = 4 . The ingredients x, y, and w are reported for different choices of the β and α parameters.
Figure 5. Electronic density (a) and KE ingredients (bd) for a jellium sphere with N = 254 electrons and r s = 4 . The ingredients x, y, and w are reported for different choices of the β and α parameters.
Computation 10 00030 g005
Figure 6. Exact Pauli Kohn–Sham enhancement factor (Equation (75)) vs. the w ingredient as computed on all grid points for several jellium spheres (see text for details). Each panel corresponds to a different choice of the β and α parameters.
Figure 6. Exact Pauli Kohn–Sham enhancement factor (Equation (75)) vs. the w ingredient as computed on all grid points for several jellium spheres (see text for details). Each panel corresponds to a different choice of the β and α parameters.
Computation 10 00030 g006
Table 1. Considered values of α and β , with the corresponding values of G 0 , see Equation (69), G p , see Equation (67), the error defined in Equation (63), and the value of 1 / F at η = 1 .
Table 1. Considered values of α and β , with the corresponding values of G 0 , see Equation (69), G p , see Equation (67), the error defined in Equation (63), and the value of 1 / F at η = 1 .
α β G 0 G p σ 1 / F ¯ ( η = 1 )
1.36111.480.07910.385
3.3114.971.480.10490.617
3.312−0.310.850.06760.545
3.312/31.862.450.07580.567
2.345/912.980.05880.469
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fabiano, E.; Sarcinella, F.; Constantin, L.A.; Della Sala, F. Kinetic Energy Density Functionals Based on a Generalized Screened Coulomb Potential: Linear Response and Future Perspectives. Computation 2022, 10, 30. https://doi.org/10.3390/computation10020030

AMA Style

Fabiano E, Sarcinella F, Constantin LA, Della Sala F. Kinetic Energy Density Functionals Based on a Generalized Screened Coulomb Potential: Linear Response and Future Perspectives. Computation. 2022; 10(2):30. https://doi.org/10.3390/computation10020030

Chicago/Turabian Style

Fabiano, Eduardo, Fulvio Sarcinella, Lucian A. Constantin, and Fabio Della Sala. 2022. "Kinetic Energy Density Functionals Based on a Generalized Screened Coulomb Potential: Linear Response and Future Perspectives" Computation 10, no. 2: 30. https://doi.org/10.3390/computation10020030

APA Style

Fabiano, E., Sarcinella, F., Constantin, L. A., & Della Sala, F. (2022). Kinetic Energy Density Functionals Based on a Generalized Screened Coulomb Potential: Linear Response and Future Perspectives. Computation, 10(2), 30. https://doi.org/10.3390/computation10020030

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop