Jellium-with-gap model applied to semilocal kinetic functionals

TL;DR

This paper introduces a new non-empirical GGA kinetic energy functional based on the jellium-with-gap model, which improves the description of weakly interacting molecular systems in DFT.

Contribution

It develops a band-gap-dependent gradient expansion and a GGA kinetic energy functional that enhances accuracy for molecular systems.

Findings

Performs well for large atoms

Accurately describes weakly interacting molecules

Based on a non-empirical, gap-dependent approach

Abstract

We investigate a highly-nonlocal generalization of the Lindhard function, given by the jellium-with-gap model. We find a band-gap-dependent gradient expansion of the kinetic energy, which performs noticeably well for large atoms. Using the static linear response theory and the simplest semilocal model for the local band gap, we derive a non-empirical generalized gradient approximation (GGA) of the kinetic energy. This GGA kinetic energy functional is remarkably accurate for the description of weakly interacting molecular systems within the subsystem formulation of Density Functional Theory.

Tables

Table 1. Table 1: Comparison of several linear-response KE gradient expansions. All integrations are performed over the volume V 𝑉 V , defined by − 1 ≤ q ≤ 1 1 𝑞 1 -1\leq q\leq 1 and Δ ≤ 1 Δ 1 \Delta\leq 1 . We show the exact KE ( T s e ​ x ​ a ​ c ​ t superscript subscript 𝑇 𝑠 𝑒 𝑥 𝑎 𝑐 𝑡 T_{s}^{exact} ) and the errors E s a ​ p ​ p ​ r ​ o ​ x = T s a ​ p ​ p ​ r ​ o ​ x − T s e ​ x ​ a ​ c ​ t superscript subscript 𝐸 𝑠 𝑎 𝑝 𝑝 𝑟 𝑜 𝑥 superscript subscript 𝑇 𝑠 𝑎 𝑝 𝑝 𝑟 𝑜 𝑥 superscript subscript 𝑇 𝑠 𝑒 𝑥 𝑎 𝑐 𝑡 E_{s}^{approx}=T_{s}^{approx}-T_{s}^{exact} (in Hartree). The GAP4 and LGAP functionals are defined in Eq. ( 17 ) and Eq. ( 20 ), respectively. The best result of each line is shown in bold style. We use LDA orbitals and densities.

atom	$T_s^{exact}$	$E_s^{GE2}$	$E_s^{Lind4}$	$E_s^{GAP4}$	$E_s^{LGAP-GE}$
Ne	125.8	-0.2	1.2	4.0	1.9
Ar	512.2	4.3	8.1	12.0	11.4
Kr	2742.3	-21.1	-6.7	1.0	10.3
Xe	7214.4	-50.7	-19.6	-8.9	23.2
Rn	21829.6	-146.5	-72.4	-56.8	48.9
Uuo	46259.6	-298.4	-162.3	-139.6	81.1
168 $e^{-}$	106907.1	-636.4	-369.8	-336.9	158.9
218 $e^{-}$	198077.5	-1065.9	-622.5	-579.7	318.9
290 $e^{-}$	389072.0	-1888.9	-1114.7	-1056.4	630.1

Table 2. Table 2: The coefficients c 0 subscript 𝑐 0 c_{0} , c 1 subscript 𝑐 1 c_{1} , and c 2 subscript 𝑐 2 c_{2} of the large- Z 𝑍 Z expansion of the kinetic energy (see Eq. ( 19 )).

	$c_0$	$c_1$	$c_2$
Exact	0.768745	-0.500	0.270
GE2	0.768745	-0.536	0.336
LGAP-GE	0.768745	-0.500	0.283

Table 3. Table 3: Mean absolute relative errors (MARE) of the non-self-consistent benchmark tests, and mean absolute errors (MAE in mHa) of FDE self-consistent tests, given by several KE functionals. The best result of each group is highlighted in bold style.

	GE2	revAPBEk	OL1	LC94	LGAP
Total KE (non-self-consistent calculations)
Atoms and ions	1.1	1.2	1.1	0.8	1.1
Jellium clusters	1.0	0.8	1.0	0.9	0.8
Jellium slabs	0.6	0.5	0.5	0.5	0.4
Molecules	0.9	0.4	0.7	0.5	0.2
KE differences (non-self-consistent calculations)
Jellium cluster DKE	27.2	23.1	28.9	21.3	22.6
Jellium surfaces	3.3	3.6	3.4	3.8	3.1
Jellium slabs dKE	5.0	3.5	4.7	4.1	3.4
Molecules AKE	184	155	185	154	159
FDE results for molecular systems 111Embedding energy errors $\Delta E=E^{\text{FDE}}-E^{\text{KS}}$ (mHa) for different KE functionals and complexes. In the last line, the mean absolute error (MAE) is reported. (self-consistent calculations)
Weak-interactions (WI)	2.46	0.13	2.49	0.36	0.21
Dipole interactions (DI)	6.48	0.48	6.59	0.67	0.45
Hydrogen bonds (HB)	10.68	1.27	10.90	1.34	1.69
Dihydrogen bonds (DHB)	4.39	3.08	4.50	2.92	2.58
Charge transfer (CT)	5.05	2.61	6.94	2.79	1.95
MAE FDE	5.66	1.72	6.31	1.82	1.50

Equations

$$T_s^{exact}=\frac{1}{2}\sum _i^{occ.}\int f_i|\nabla {\varphi }_i(\mathbf{r}){|}^2{\mathrm{d}}^3\mathbf{r}.$$

$$F^{Lind}={\left(\frac{1}{2}+\frac{1-{\eta }^2}{4\eta }\ln \left|\frac{1+\eta }{1-\eta }\right|\right)}^{-1},$$

$$-\frac{1}{{\chi }^{Jell.}}=\frac{{\pi }^2}{k_F}{F}^{Lind}.$$

$$-\frac{1}{{\chi }^{Semic.}(k)}\underset{k\to 0}{⟶}\frac{b}{k^2},$$

$$F^{Lind}$$

$$T_s^{Lind4}[n]=\int d\mathbf{r}{\tau }^{TF}\left(1+\frac{5}{27}{s}^2+\frac{8}{81}{q}^2\right),$$

$$T_s^{GE2}[n]=\int d\mathbf{r}{\tau }^{TF}(1+\frac{5}{27}{s}^2)$$

$$T_s^{GE4}[n]=\int d\mathbf{r}{\tau }^{TF}(1+\frac{5}{27}{s}^2+\frac{8}{81}{q}^2-\frac{1}{9}{s}^{2}q+\frac{8}{243}{s}^4),$$

$$1/F^{GAP}=\frac{1}{2}-\frac{\Delta (\arctan (\frac{4\eta +4{\eta }^2}{\Delta })+\arctan (\frac{4\eta -4{\eta }^2}{\Delta }))}{8\eta }+$$

$$+(\frac{{\Delta }^2}{128{\eta }^3}+\frac{1}{8\eta }-\frac{\eta }{8})\ln (\frac{{\Delta }^2+(4\eta +4{\eta }^2{)}^2}{{\Delta }^2+(4\eta -4{\eta }^2{)}^2}),$$

$$F^{GAP}⟶\frac{3{\Delta }^2}{16{\eta }^2}+\frac{9}{5}+\frac{3}{175}\frac{175{\Delta }^2-192}{{\Delta }^2}{\eta }^2+$$

$$-\frac{64}{875}\frac{525{\Delta }^2-368}{{\Delta }^4}{\eta }^4+\mathcal{O}({\eta }^6)\mathrm{when}\eta \to 0.$$

$$F^{GAP}⟶\left[1+\frac{1}{3}{\eta }^2+\frac{8}{45}{\eta }^4+...\right]+$$

$$\displaystyle\Delta\left[\frac{\pi}{8}\frac{1}{\eta}+\frac{\pi}{12}\eta+\frac{7\pi}{120}\eta^{3}+...\right]+\Delta^{2}\Bigg{[}\frac{\pi^{2}-4}{64}\frac{1}{\eta^{2}}+\frac{3\pi^{2}-16}{192}+$$

$$\displaystyle\left(\frac{-17}{180}+\frac{13\pi^{2}}{960}\right)\eta^{2}+\left(\frac{-383}{3780}+\frac{683\pi^{2}}{60480}\right)\eta^{4}+...\Bigg{]}+...\;.$$

$$F^{GAP}=F^{Lind},\mathrm{when}\mathrm{\Delta }=0.$$

$$\underset{\Delta \to 0}{\lim }\underset{\eta \to 0}{\lim }{F}^{GAP}$$

$$\underset{\eta \to 0}{\lim }\underset{\Delta \to 0}{\lim }{F}^{GAP}$$

$$F^{GAP}\to 3{\eta }^2-\frac{3}{5}+$$

$$(-\frac{24}{175}+\frac{3}{16}{\Delta }^2)\frac{1}{{\eta }^2}+\mathcal{O}(\frac{1}{{\eta }^4}).$$

$$T_s^{GAP4}[n]=\int d\mathbf{r}{\tau }^{TF}(\frac{a_1}{s^2}+\frac{a_2}{s}+a_3+a_{4}s+a_{5}q+$$

$$+a_6{s}^2+a_{7}sq+a_8{s}^3+a_9{s}^4+a_{10}{q}^2+a_{11}{s}^{2}q).$$

$$F(\eta )=\frac{k_F}{{\pi }^2}\mathcal{F}\left(\frac{{\delta }^2{T}_s[n]}{\delta n(\mathbf{r})\delta n({\mathbf{r}}^{\prime })}{|}_{n_0}\right),$$

$$T_s^{GAP4}[n]=\int d\mathbf{r}{\tau }^{TF}[{\Delta }^2\frac{27}{91}\frac{{\pi }^2-4}{64}\frac{1}{s^2}+\Delta \frac{5\pi }{72}\frac{1}{s}+1+$$

$${\Delta }^2(\frac{{\pi }^2}{64}-\frac{1}{12})+\Delta \frac{5\pi }{36}s+(\frac{5}{27}+{\Delta }^2(\frac{-17}{324}+\frac{13{\pi }^2}{1728}))s^2+$$

$$\Delta \frac{-7\pi }{216}sq+(\frac{8}{81}+{\Delta }^2(\frac{-383}{6804}+\frac{683{\pi }^2}{108864}))q^2].$$

$$E_g(\mathbf{r})=a|\nabla n(\mathbf{r}){|}^m/n(\mathbf{r}{)}^{2(2m-1)/3},m\ge 0,a\ge 0.$$

$$T_s=c_0{Z}^{7/3}+c_1{Z}^2+c_2{Z}^{5/3}+...,$$

$$T_s^{LGAP-GE}=\int d\mathbf{r}{\tau }^{TF}[1+a\frac{5\pi }{72}s+(\frac{5}{27}+a^2\frac{27}{91}\frac{{\pi }^2-4}{64})s^2+$$

$$a\frac{5\pi }{36}{s}^3+\mathcal{O}(|\nabla n{|}^4)]$$

$$=\int d\mathbf{r}{\tau }^{TF}[1+0.0131s+0.18528s^2+0.0262s^3].$$

$$F_s^{OL1}=1+\frac{5}{27}{s}^2+cs,$$

$$\frac{\delta T_s}{\delta n}=\frac{\partial {\tau }^{TF}}{\partial n}{F}_s(s)+{\tau }^{TF}\frac{\partial F_s}{\partial s}\frac{\partial s}{\partial n}-\nabla \cdot [\frac{1}{s}\frac{\partial F_s}{\partial s}\cdot \frac{\nabla n}{n^{8/3}}],$$

$$F^{LGAP}_{s}=1+\kappa\Big{(}1-e^{-\mu_{1}s-\mu_{2}s^{2}-\mu_{3}s^{3}}\Big{)},$$

Full text

I Introduction

Density Functional Theory (DFT) Parr and Yang (1989); Dreizler and Gross (1990) is the most used computational method for electronic structure calculations of molecular and extended systems, providing the highest accuracy/computational cost ratio. In the conventional DFT formalism, the Kohn-Sham (KS) scheme Kohn and Sham (1965), the ground-state electronic density $n(\mathbf{r})$ is determined from a set of auxiliary KS orbitals ($\phi_{i}(\mathbf{r})$): the KS-DFT method is exact but for the approximations of the exchange-correlation (XC) functional. However, for large scale calculations, the computational cost of KS-DFT becomes unaffordable, as one needs to compute all the occupied KS orbitals in order to construct the density as $n(\mathbf{r})=\sum_{i}^{occ.}f_{i}|\phi_{i}(\mathbf{r})|^{2}$, where $f_{i}$ is the occupation number (2, for closed-shell systems).

Among other linear scaling methods Yang (1991); Goedecker (1999); Bowler et al. (2002); Kussmann et al. (2013), two DFT methods are attracting strong interest: i) In the orbital-free version of DFT (OF-DFT) Smargiassi and Madden (1994); Wang and Carter (2002); Gavini et al. (2007); Wesolowsky and Wang (2013); Chen et al. (2016), $n(\mathbf{r})$ can be computed directly via the Euler equation Parr and Yang (1989), without the need of KS orbtials; ii) In the subsystem version of DFT (Sub-DFT) Cortona (1991); Wesolowski (2006); Krishtal et al. (2015); Wesolowski et al. (2015), also known as Frozen-Density-Embedding (FDE), $n(\mathbf{r})$ is computed as the sum of the electronic densities of several (smaller) subsystems in which the total system is partitioned, which can be computed simultaneously. Both approaches allow in principle calculations of large systems, but the final accuracy depends directly on the approximations of the non-interacting kinetic energy (KE) functional $T_{s}$ (which are definitely more important than the ones for the XC energy, that are also present in standard KS calculations). We recall that the exact KS KE functional is:

[TABLE]

Thus the KE is explicitly known only as a function of $\phi_{i}$ but not as a functional of $n$.

On the other hand, in Sub-DFT the interaction between the subsystems is taken into account via the so called embedding potentials Wesolowski (2006); Krishtal et al. (2015); Wesolowski et al. (2015), which depends on the non-additive-KE: in the case of just two subsystems (A and B) it is $T_{s}^{nadd}[n_{A};n_{B}]=T_{s}[n_{A}+n_{B}]-T_{s}[n_{A}]-T_{s}[n_{B}]$.

The development of an accurate approximation of $T_{s}[n]$ (and/or $T_{s}^{nadd}[n_{A};n_{B}]$) is one of the biggest DFT challenges Karasiev and Trickey (2012); Xia et al. (2012); Lastra et al. (2008). Nowadays, the most sophisticated KE approximations have been constructed to be exact for the linear response of jellium model, by incorporating the Lindhard function in their fully non-local expressions Wang and Carter (2002); Wang et al. (1998); Ho et al. (2008); García-González et al. (1996, 1998); Garcia-Aldea and Alvarellos (2008). We recall that the Lindhard function Lindhard (1954); Wang and Carter (2002)

[TABLE]

where $\eta=k/(2k_{F})$ is the dimensionless momentum ($k_{F}=(3\pi^{2}n)^{1/3}$ being the Fermi wave vector of the jellium model with the constant density $n$), is related to the Jellium density response $\chi^{Jell}$ viaWang and Carter (2002)

[TABLE]

The non-local KE functionals based on the Lindhard function are accurate for simple metals where the nearly-free electron gas is an excellent model but they can not describe well semiconductors and insulators, where the density response function behaves as Pick et al. (1970); Huang and Carter (2010)

[TABLE]

with $b$ being positive and material-dependent. Several KE functionals have been constructed to improve the description of semiconductors Huang and Carter (2010); Shin and Carter (2014), but Eq. (4) has not been explicitly used in their expressions due to the lack of a sophisticated analytical form that can recover both the Lindhard function and Eq. (4).

In this article, we will investigate the generalization of the Lindhard function for the jellium-with-gap model which satisfies Eq. (4).

The jellium-with-gap modelRey and Savin (1998), was developed outside the KS framework, using perturbation theory to take into account the band gap energy. This model was used to have qualitative and quantitative insight for semiconductors Callaway (1959); Penn (1962); Srinivasan (1969); Levine and Louie (1982); Tsolakidis et al. (2004), to develop an XC kernel for the optical properties of materials Trevisanutto et al. (2013), and to construct accurate correlation energy functionals for the ground-state DFT Rey and Savin (1998); Krieger et al. (1999, 2001); Toulouse et al. (2002); Toulouse and Adamo (2002); Fabiano et al. (2014). We will show that the Lindhard function for the jellium-with-gap model ($F^{GAP}$), previously introduced by Levine and Louie Levine and Louie (1982) in a different context (dielectric constant and XC potential), may be seen as a sophisticated analytical form suitable for KE approximations.

The article is organized as follow:

In Section II we discuss the properties of $F^{GAP}$, we derive its (band-gap-dependent) KE gradient expansion, and we assess it for large atoms. By using a local gap model, we propose a simple KE gradient expansion that it is very accurate for the semiclassical atom theory.

In Section III we discuss the implications of this result in DFT by constructing a simple KE functional at the Generalized Gradient Approximation (GGA) level of theory based on the gradient expansion of the jellium-with-gap model. GGA KE functionals are computationally very efficient and play a key role for the simulation of large systems. We mention that the development of semilocal KE functionals is nowadays an active field Constantin et al. (2011); Laricchia et al. (2011a); Karasiev et al. (2013); Laricchia et al. (2013a); Borgoo and Tozer (2013); Borgoo et al. (2014); Lindmaa et al. (2014); Della Sala et al. (2015); Cancio et al. (2016); Śmiga et al. (2017).

Finally, in Section IV we summarize our results.

II Theory

II.1 Properties and gradient expansions for the jellium model

For the conventional infinite jellium model, the Lindhard function behaves as:

[TABLE]

Equation (5) contains important physics that has been used in the construction of semilocal KE density functionals Wang and Carter (2002). Thus, the KE gradient expansion which recovers the first three terms in the right hand side of Eq. (5) can be easily derived Wang and Carter (2002) (see also Eqs. (15) and (16) in section II-C and the corresponding discussion). It is

[TABLE]

where $\tau^{TF}=\frac{3}{10}(3\pi^{2})^{2/3}n^{5/3}$ is the Thomas-Fermi KE density Thomas (1927); Fermi (1927), which is exact for the jellium model, and $s=|\nabla n|/[2k_{F}n]$, $q=\nabla n^{2}/[4(3\pi^{2})^{2/3}n^{5/3}]$ are the reduced gradient and Laplacian, respectively. Equation (7) resembles the second-order gradient expansion Kirzhnitz (1957) (GE2)

[TABLE]

(derived also within the linear response of the jellium model), as well as the fourth-order gradient expansion Hodges (1973); Brack et al. (1976); Laricchia et al. (2013a) of the KE

[TABLE]

with the exception of the terms $\propto s^{2}q$, $\propto s^{4}$, which are beyond the linear response.

Note that $F^{Lind}(\eta=0)=1$ is the leading term in the expansion of Eq. (5) and it corresponds to the Thomas-Fermi local density approximation, whose linear response in the wave vector space is just the Fourier transform of the second-functional derivative, i.e. $\delta^{2}T_{s}^{TF}/\delta n(\mathbf{r})\delta n(\mathbf{r}^{\prime})\sim k_{F}^{-1}\delta(\mathbf{r}-\mathbf{r}^{\prime})$. We recall that the limit $\eta=0$ is very powerful, being also used in the construction of the adiabatic local density approximation (ALDA) XC kernel of the linear response time-dependent DFT Burke et al. (2005); Constantin and Pitarke (2007).

II.2 Properties of the Lindhard function for the jellium-with-gap model

Levine and Louie Levine and Louie (1982) proposed the density-response function $\chi^{GAP}(k,\omega)$ of the jellium-with-gap model, and the corresponding [i.e. from Eq. (3)] Lindhard function for jellium-with-gap model is

[TABLE]

where $\Delta=2E_{g}/k_{F}^{2}$ and $E_{g}$ is the gap.

For a given $\Delta$, a series expansion of $F^{GAP}$ for $\eta\rightarrow 0$ gives:

[TABLE]

Thus, for any system with $\Delta>0$ we have that $F^{GAP}\propto\Delta^{2}\eta^{-2}$. This term is correct (see Eq. (4)) and it has been also used in the jellium-with-gap XC kernel Trevisanutto et al. (2013), which gives accurate optical absorption spectra of semiconductors and insulators. On the other hand, if we first perform a series expansion for $\Delta\rightarrow 0$, and then a series expansion for $\eta\rightarrow 0$ we obtain:

[TABLE]

Equation (10) confirms that, by construction, we have

[TABLE]

Inspection of Eqs. (9) and (10) clearly shows that

[TABLE]

meaning that $F^{GAP}$ has an “order of limits problem”. Such a situation is common in DFT. For example, we recall that several meta-GGA XC functionals (e.g. TPSS Tao et al. (2003), revTPSS Perdew et al. (2009, 2011), BLOC Constantin et al. (2013a, b), SA-TPSS Constantin et al. (2016a), VT{8,4} del Campo et al. (2012)) suffer of such a order of limits problem. Nonetheless, they are accurate for many systems and properties, showing realistic system-averaged XC hole models Constantin et al. (2013b).

In the opposite limit, i.e. for $\eta\rightarrow\infty$, we have

[TABLE]

Therefore, in this limit, $F^{GAP}$ always behaves as $F^{Lind}$ for $\Delta=0$.

In the upper panel of Fig. 1, we show $1/F^{GAP}$ for several values of $\Delta$. The plots are all smooth. At large $\eta$, $F^{GAP}$ recovers the Lindhard function [see Eq. (14)], while at small $\eta$ it is driven by the term $\propto\eta^{-2}$. The plot of the linear response of $T_{s}^{Lind4}[n]$ (Eq. (5)) is also given for comparison.

In the lower panel of Fig. 1 we report the accuracy of Eq. (10), considering only the terms explicitly indicated in the equation, for $\Delta=0.1,0.5,$ and $1$. Even for the case $\Delta=1$, this expansion is still very accurate for $\eta\leq 1$.

II.3 Kinetic energy gradient expansions from the linear-response of the jellium-with-gap model

Next we proceed to build the linear-response jellium-with-gap KE gradient expansion, that should recover Eq. (7) when $\Delta=0$. To this purpose, we consider the GAP4 expansion, with the general form of the KE fourth-order gradient expansion

[TABLE]

Performing the linear response of such a functional

[TABLE]

where $\mathcal{F}$ represents the Fourier transform, we can find the coefficients $a_{i}$, by comparing term-by-term with Eq. (10). Nevertheless, the straightforward calculation of Eq. (15) requires a tedious and long algebra Nazarov and Vignale (2011); Karimi and Ullrich (2014). Instead, a more elegant and simpler way to obtain the linear response of a given semilocal functional has been proposed in Ref. Tao et al., 2008: consider a small pertubation in density at $\mathbf{r}=\mathbf{0}$, of the form $n=n_{0}+n_{k}e^{i\mathbf{k}\mathbf{r}}$, such that $\nabla n=n_{k}i\mathbf{k}e^{i\mathbf{k}\mathbf{r}}$, and $\nabla^{2}n=-n_{k}k^{2}e^{i\mathbf{k}\mathbf{r}}$, with $n_{k}\ll n_{0}$. Thus, at $\mathbf{r}=\mathbf{0}$, these expressions are simply $n=n_{0}+n_{k}$, $\nabla n=n_{k}i\mathbf{k}$, and $\nabla^{2}n=-n_{k}k^{2}$. Inserting them in the functional expression, the linear response is obtained as twice the second-order coefficient of the series expansion with respect to $n_{k}/n_{0}$. After some algebra, the KE gradient expansion which gives the linear response of Eq. (10) is found to be

[TABLE]

The terms $\propto s^{-2}$ and $\propto s^{-1}$ account for the terms $\propto\eta^{-2}$ and $\propto\eta^{-1}$ of Eq. (10). These terms contribute only for a non-zero gap, i.e. in semiconductors and insulators, but not in metals. At $\Delta=0$, $T_{s}^{GAP4}[n]$ correctly recovers $T_{s}^{Lind4}[n]$.

To test $T_{s}^{GAP4}[n]$, we perform calculations for noble atoms, up to $Z=290$ electrons, using LDA orbitals and densities, in the Engel code Engel and Vosko (1993); Engel (2003). We consider $\Delta=2E_{g}/k_{F}^{2}(\mathbf{r})$ with $E_{g}$ being the KS band gap of the atoms. Because the gradient expansion is well defined only at small gradients and small-$\Delta$, we perform all the integrations over the volume $V$ defined by the conditions $-1\leq q\leq 1$ and $\Delta\leq 1$, in a similar manner as in Ref. Constantin et al., 2016b.

The results are reported in Table 1. For small atoms (Ne and Ar), the GE2 is more accurate than $T_{s}^{Lind4}[n]$ and $T_{s}^{GAP4}[n]$. However, we recall that in the case of a small number of electrons, the semiclassical and statistical concepts beyond the gradient expansions do not hold. In fact, for larger atoms (Kr to the noble atom with 290 electrons), both $T_{s}^{GAP4}$ and $T_{s}^{Lind4}$ outperform GE2. In particular, $T_{s}^{GAP4}$ shows the best performance, improving over $T_{s}^{Lind4}[n]$ and proving that, due to the inclusion of the gap, $F^{GAP4}$ contains important physics beyond $F^{Lind4}$.

II.4 Local band-gap

In order to use Eq. (17) in semilocal DFT, we need to replace the true band gap $E_{g}$, with a density dependent local band gap. There are several models for the local band gap Fabiano et al. (2014); Krieger et al. (1999), constructed from the exponentially decaying density behavior Krieger et al. (1999) or from conditions of the correlation energy Fabiano et al. (2014). In the slowly-varying density limit, they behave as $E_{g}\sim|\nabla n|^{m}$, with $m\geq 2$. However, none of them can be considered accurate in this density regime.

On the other hand, under a uniform density scaling $n_{\lambda}(\mathbf{r})=\lambda^{3}n(\lambda\mathbf{r})$, the local band gap should behave as $E_{g}\sim\lambda^{2}$. This condition is fulfilled by the general formula

[TABLE]

Because other exact conditions of the local gap in the slowly varying density limit are not known, we use Eq. (18) in the expression of $T_{s}^{GAP4}$, considering the case with $m=2$. We fix the parameter $a$ requiring that the gradient expansion should recover the first two terms of the kinetic energy asymptotic expansion for the large, neutral atom Thomas (1927); Fermi (1927); Scott (1952); Schwinger (1980, 1981); Englert and Schwinger (1984, 1985); Elliott et al. (2008); Lee et al. (2009); Fabiano and Constantin (2013)

[TABLE]

where $Z$ is the number of electrons. The first term in Eq. (19) is the Thomas-Fermi one Thomas (1927); Fermi (1927), the second is the Scott correction due to the atomic inner core Scott (1952), and the last term accounts for quantum oscillations Schwinger (1980, 1981); Englert and Schwinger (1984, 1985). The exact coefficients are shown in the first line of Table 2. As in Ref. Lee et al., 2009, we assume that any gradient expansion that is exact for the uniform electron gas, should have the exact $c_{0}$ coefficient. The calculation of $c_{1}$ and $c_{2}$ has been done using the method proposed in Ref. Lee et al., 2009. We recall that the semiclassical atom theory has been often used in the development of exchange functionals Constantin et al. (2016b, 2011); Sun et al. (2015a, b); Becke (1988); Cancio et al. (2012) and occasionally also for kinetic energy functionals Laricchia et al. (2011a). Finally, we mention that these gradient expansions are models for the total KE, and not for the KE density, where the use of the reduced Laplacian $q$ (which is not present in linear response of the jellium model) is essential Cancio et al. (2016); Cancio and Redd (2016); Lindmaa et al. (2014); Śmiga et al. (2015).

Using the procedure described above, we find $a=0.0075$, and we obtain the following gradient expansion (denoted as LGAP-GE)

[TABLE]

Note that in Eq. (20) only terms to up $s^{3}$ are considered (terms in Eq. (17) proportional to $q$ or $q^{2}$ are neglected, as these terms will correspond to $s^{4}$).

As shown in Table 2, LGAP-GE gives a very accurate large-$Z$ expansion, having the $c_{2}$ coefficient close to exact. The results for noble atoms are reported in Table 1. LGAP-GE is reasonably accurate for all atoms and, as expected due to the inclusion of the semiclassical atom theory, the accuracy increases with the number of electrons.

One additional observation is that LGAP-GE contains odd powers of the reduced gradient, in contrast with $F^{Lind4}$. Nevertheless, Ou-Yang and Mel Levy have already shown that using non-uniform coordinates scaling requirements HuiOu-Yang and Levy (1990), the GE4 terms in the KE gradient expansion can be replaced by an $s$-only dependent term Ou-Yang and Levy (1991), whose coefficient must be positive (and was fitted to the Xe atom). The resulting simple KE functional, that behaves better than GE4 for the non-uniform density scaling, has the following enhancement factor ($F_{s}=\tau^{approx}/\tau^{TF}$):

[TABLE]

with $c=0.01459$ being slightly bigger than its LGAP-GE counterpart. Anyway, we need to acknowledge that, since the kinetic potential of a GGA functional (with the enhancement factor $F_{s}$) has the general form

[TABLE]

a necessary condition for it to be well defined is $|\frac{1}{s}\frac{\partial F_{s}}{\partial s}|<\infty$. This is not satisfied by the LGAP-GE (and OL1). Thus, the term $\propto s$ gives a diverging kinetic potential ($\delta T_{s}/\delta n\rightarrow\infty$) at $s=0$. This is due to the high non-locality of Eq. (17), which was not fully suppressed by the local gap model of Eq. (18) with $m=2$. Note that this divergence is a direct consequence of the jellium-with-gap theory. Nevertheless, for molecular systems $s=0$ only at the middle of bonds, and it has been found that this divergence is not important in real calculations of weakly-bounded molecular systems Gotz et al. (2009). In fact, the same problem is shared by other well-known KE functionals Ou-Yang and Levy (1991); Thakkar (1992); Borgoo and Tozer (2013).

III Kinetic energy functional constructed from the LGAP gradient expansion

III.1 The LGAP GGA

To show the importance of the LGAP-GE, we construct a simple GGA functional (named LGAP-GGA or simply LGAP) that recovers the LGAP-GE in the slowly-varying density regime. We consider the RPBE exchange enhancement factor form Hammer et al. (1999), $F_{x}^{RPBE}=1+\kappa(1-e^{-\mu s^{2}/\kappa})$, and we fix $\kappa=0.8$ from the Lieb-Oxford bound Lieb and Oxford (1981), using the approximate link between the kinetic and exchange energies (i.e. the conjointness conjecture Lee et al. (1991); March and Santamaria (1991); Constantin et al. (2011)). Note that, to our knowledge, the RPBE functional form has not been yet used in the development of kinetic functionals. The LGAP kinetic enhancement factor is therefore defined as

[TABLE]

where $\mu_{1}=b_{1}/\kappa$, $\mu_{2}=b_{2}/\kappa+\mu_{1}^{2}/2$, and $\mu_{3}=b_{3}/\kappa+\mu_{1}\mu_{2}-\mu_{1}^{3}/6$, such that it recovers the LGAP-GE in the slowly-varying density limit. Here $b_{1}=0.0131$, $b_{2}=0.18528$, and $b_{3}=0.0262$ [see Eq. (20)].

III.2 The kinetic energy benchmark

In order to assess the LGAP KE functional, we consider several known tests.

For total KE:

• •

The benchmark set of atoms and ions Perdew and Constantin (2007); Laricchia et al. (2013a, 2011a). All calculations employed analytic Hartree-Fock orbitals and densities Clementi and Roetti (1974);

• •

The Na jellium clusters ($r_{s}=3.93$) set with magic electron numbers 2, 8, 18, 20, 34, 40, 58, 92, and 106, used in Refs. Perdew and Constantin, 2007; Laricchia et al., 2013a, 2011a. We use exact exchange orbitals and densities;

• •

The set of two interacting jellium slabs at different distances Laricchia et al. (2013a). Each jellium slab has $r_{s}=3$ and a thickness of $2\lambda_{F}$. Here $\lambda_{F}=2\pi/k_{F}$ is the Fermi wavelength. The calculations were performed using the orbitals and densities resulting from numerical Kohn-Sham calculations within the local density approximation Kohn and Sham (1965) for the XC functional;

• •

The set of molecules (H2, HF, H2O, CH4, NH3, CO, F2, HCN, N2, CN, NO, and O2) used in Refs. Iyengar et al., 2001; Perdew and Constantin, 2007; Laricchia et al., 2013a. The noninteracting kinetic energies of test molecules were calculated using the PROAIMV code Biegler-könig et al. (1982). The required Kohn-Sham orbitals were obtained by Kohn-Sham calculations performed with the uncontracted 6-311+G(3df,2p) basis set, the Becke 1988 exchange functional Becke (1988), and Perdew-Wang correlation functional Perdew and Wang (1992).

For KE differences:

• •

The disintegration kinetic energy (DKE) of a jellium cluster Constantin and Ruzsinszky (2009); Laricchia et al. (2013a);

• •

The jellium surfaces test with bulk parameter $r_{s}$=2, 4, and 6 into the liquid drop model (LDM), as in Refs. Perdew and Constantin, 2007; Laricchia et al., 2013a, 2011a;

• •

The dissociation KE (dKE) of a jellium slab into two pieces (as in Ref. Laricchia et al., 2013a);

• •

The atomization KE (AKE) of molecules Iyengar et al. (2001); Perdew and Constantin (2007); Laricchia et al. (2013a).

For non-additive KE:

We employ the LGAP functional in subsystem DFT calculations, using the TURBOMOLE tur program, together with FDE script Laricchia et al. (2010). The FDE calculations have been performed with a supermolecular def2-TZVPPD Weigend and Ahlrichs (2005); Rappoport and Furche (2010) basis set and the Perdew-Burke-ErnzerhofPerdew et al. (1996) XC functional. Five weakly interacting groups of molecular complexes are considered as a benchmark Laricchia et al. (2011a, 2013a, b, 2012, 2013b); Śmiga et al. (2015):

• WI:

weak interaction (He-Ne, He-Ar, Ne2, Ne-Ar, CH4-Ne, C6H6-Ne, (CH4)2);

• DI:

dipole-dipole interaction ((H2S)2, (HCl)2, H2S-HCl, CH3Cl-HCl,CH3SH-HCN, CH3SH-HCl);

• HB:

hydrogen bond ((NH3)2, (HF)2, (H2O)2, HF-HCN, (HCONH2)2, (HCOOH)2);

• DHB:

double hydrogen bond (AlH-HCl, AlH-HF, LiH-HCl, LiH-HF, MgH2-HCl, MgH2-HF, BeH2-HCl, BeH2-HF);

• CT:

charge transfer (NF3-HCN,C2H4-F2,NF3-HCN, C2H4-Cl2, NH3-F2, NH3-ClF, NF3-HF, C2H2-ClF, HCN-ClF, NH3-Cl2, H2O-ClF, NH3-ClF).

III.3 Results

We compare our results with revAPBEk Constantin et al. (2011) and LC94 Lembarki and Chermette (1994) GGAs, which are considered state-of-the-art KE functionals for FDE Laricchia et al. (2011a), as well as with GE2 Kirzhnitz (1957); Brack et al. (1976) and OL1 Ou-Yang and Levy (1991). The KE enhancement factors of the considered functionals are reported in Fig. 2. In the inset of Fig. 2, we show that LGAP and LGAP-GE are identical (by construction) at relatively small values of the reduced gradient ($0\leq s\leq 0.5$), both differing significantly from the GE2 behavior. Consequently, LGAP shows a bigger enhancement factor than both LC94 and revAPBEk (i.e. $F_{s}^{LGAP}\geq\sim F_{s}^{revAPBEk}\geq\sim F_{s}^{LC94}$) when $s\leq 2.5$. Such a feature has been proved to be essential for jellium surfaces Vitos et al. (2000). On the other hand, the LGAP enhancement factor recovers its maximum value $F_{s}\rightarrow 1+\kappa$ at $s\approx 3$, faster than revAPBEk.

In Table 3 we report the numerical results of all the tests. For total KE tests, LGAP gives the best overall performance, among the considered functionals, being the best for jellium clusters, jellium slabs and molecules. For KE differences LGAP is the most accurate for jellium surfaces and dissociation KE of jellium slabs. We also mention that LGAP performs reasonably well for all the other tests, being in line with revAPBEk.

Finally, LGAP outperforms the other functionals for the FDE theory, being especially accurate for dipole-dipole, dihydrogen bond and charge transfer interactions. These latter results show that, in agreement with the finding of Ref. Gotz et al., 2009, the divergence at $s=0$ of the LGAP kinetic potential is not important for calculations of weakly-bounded molecular systems. Moreover, results indicate that the LGAP-GE gradient expansion can be successfully used in the kinetic energy functional construction, which perform relatively well in FDE theory.

IV Conclusions

In conclusion, we have investigated the linear response of the jellium-with-gap model, in the context of semilocal kinetic functionals. We have shown that the Levine and Louie Levine and Louie (1982) analytical generalization of the Lindhard function ($F^{GAP}$) contains important physics beyond jellium model, and in particular we mention the following properties:

$(i)$ $F^{GAP}$ recovers the Lindhard function when the band gap is zero (i.e. $E_{g}=0$);

$(ii)$ $F^{GAP}$ has the correct behavior (see Eq. (4)) at small wave vectors, expressing the material-dependent constant $b$ in terms of the band gap;

$(iii)$ In the regime of small band gap energy (i.e. $E_{g}\leq E_{F}$, with $E_{F}$ being the Fermi energy), $F^{GAP}$ gives the GAP4 gradient expansion of the kinetic energy (see Eq. (17)), which is band-gap-dependent, and performs remarkably well in the atomic regions where the density varies slowly, improving over $T_{s}^{Lind4}$ of Eq. (7) (see Table I).

These features show that $F^{GAP}$ should be further investigated and exploited in the field of non-local kinetic functionals Wang and Carter (2002); Huang and Carter (2010); Wang et al. (1998); Shin and Carter (2014); Ho et al. (2008); Alonso and Girifalco (1978); García-González et al. (1996); Chacón et al. (1985); García-González et al. (1998); Garcia-Aldea and Alvarellos (2008); Karasiev et al. (2014a, b); Karasiev and Trickey (2015); Cangi et al. (2010); Ribeiro and Burke (2016), and we will like to address this important issue in further work.

Finally, by considering a local band gap model, and a simple enhancement factor form, we have constructed the non-empirical LGAP GGA kinetic energy functional, derived from the linear response of the jellium-with-gap model (a.i. the GAP4 gradient expansion). This functional showed the best performance in the context of FDE theory. Thus, it can be further used in real applications.

Acknowledgements

This work was partially supported by the National Science Center under Grant No. DEC-2013/11/B/ST4/00771 and DEC-2016/21/D/ST4/00903.