Nonlocal Kinetic Energy Functionals By Functional Integration

TL;DR

This paper introduces a new nonlocal kinetic energy functional, MGP, constructed via functional integration, which accurately models electron densities and energies in various materials with low computational cost.

Contribution

The work presents a novel nonlocal kinetic energy functional derived from the second derivative, satisfying key physical conditions and enabling efficient orbital-free DFT calculations.

Findings

MGP accurately predicts equilibrium properties of silicon and III-V semiconductors.

MGP is numerically stable and converges within 12 iterations.

Computational cost of MGP is comparable to existing GGA functionals.

Abstract

Since the seminal works of Thomas and Fermi, researchers in the Density-Functional Theory (DFT) community are searching for accurate electron density functionals. Arguably, the toughest functional to approximate is the noninteracting Kinetic Energy, $T_s[\rho]$, the subject of this work. The typical paradigm is to first approximate the energy functional, and then take its functional derivative, $\frac{\delta T_s[\rho]}{\delta \rho(r)}$, yielding a potential that can be used in orbital-free DFT, or subsystem DFT simulations. Here, this paradigm is challenged by constructing the potential from the second-functional derivative via functional integration. A new nonlocal functional for $T_s[\rho]$ is prescribed (which we dub MGP) having a density independent kernel. MGP is constructed to satisfy three exact conditions: (1) a nonzero "Kinetic electron" arising from a nonzero exchange hole; (2) the second functional derivative must reduce to the inverse Lindhard function in the limit of homogenous densities; (3) the potential derives from functional integration of the second functional derivative. Pilot calculations show that MGP is capable of reproducing accurate equilibrium volumes, bulk moduli, total energy, and electron densities for metallic (BCC, FCC) and semiconducting (CD) phases of Silicon as well as of III-V semiconductors. MGP functional is found to be numerically stable typically reaching selfconsistency within 12 iteration of a truncated Newton minimization algorithm. MGP's computational cost and memory requirements are low and comparable to the Wang-Teter (WT) nonlocal functional or any GGA functional.

Tables

Table 1. Table 1: Equilibrium volume ( V 0 subscript 𝑉 0 V_{0} in Bohr 3 ), minimum energy ( E 0 subscript 𝐸 0 E_{0} in eV/atom), and bulk modulus ( B 0 subscript 𝐵 0 B_{0} in GPa) for metallic phases (BCC,FCC) and semiconductor phase (CD) of silicon computed by OF-DFT/BLPS and OF-DFT/OEPP using the MGP, WGC, WT, and HC KEDFs with optimal choices of their parameters. The corresponding benchmark KS-DFT/BLPS and KS-DFT/OEPP results are also given. HC KEDF values are taken from Ref. 10 . WGC’s optimal parameters are taken from Ref. 22 and 15 .

Systems	Functionals	Pseudopotentials	Parameters	$V_0(a.u/Cell)$	$E_0\left(eV\right)$	$B_0\left(GPa\right)$
FCC	KS	OEPP	– –	118.2	-107.761	71
	MGP	OEPP	a=0.5,b=0.3	119.5	-107.765	73
	WGC	OEPP	$\gamma =2.2,\rho \ast ={\rho }_0$	112.7	-107.733	113
	WGC	OEPP	$\gamma =1.0,\rho \ast =1.05{\rho }_0$	115.9	-107.790	98
	KS	BLPS	– –	97.0	-109.248	83
	MGP	BLPS	a=0.6,b=0.4	97.6	-109.243	75
	WT	BLPS	– –	97.5	-109.260	58
	WGC	BLPS	$\gamma =2.2,\rho \ast ={\rho }_0$	94.4	-109.258	96
BCC	KS	OEPP	– –	234.7	-215.490	74
	MGP	OEPP	a=0.66,b=0.38	235.7	-215.492	89
	WGC	OEPP	$\gamma =2.2,\rho \ast =1.05{\rho }_0$	228.9	-215.504	106
	WGC	OEPP	$\gamma =1.0,\rho \ast =1.05{\rho }_0$	230.5	-215.521	101
	KS	BLPS	– –	197.1	-218.556	98
	MGP	BLPS	a=0.94,b=0.6	205.4	-218.553	97
	WT	BLPS	– –	204.3	-218.666	66
	WGC	BLPS	$\gamma =2.2,\rho \ast =1.05{\rho }_0$	196.9	-218.553	106
CD	KS	OEPP	– –	309.3	-215.538	61
	MGP	OEPP	a=0.341,b=0.45	310.0	-215.534	54
	WGC	OEPP	$\gamma =4.2,\rho \ast =1.05{\rho }_0$	352.2	-215.512	34
	WGC	OEPP	$\gamma =5.0,\rho \ast =1.05{\rho }_0$	332.7	-214.974	44
	KS	BLPS	– –	266.9	-219.258	98
	MGP	BLPS	a=0.364,b=0.57	265.6	-219.258	95
	HC	BLPS	$\lambda =0.01,\beta =0.65$	269.4	-219.248	97
	WT	BLPS	– –	–	–	–
	WGC	BLPS	$\gamma =4.2,\rho \ast =1.05{\rho }_0$	287.4	-218.838	64

Table 2. Table 2: Equilibrium volume ( V 0 subscript 𝑉 0 V_{0} in Bohr 3 ), minimum energy ( E 0 subscript 𝐸 0 E_{0} in eV/atom), and bulk modulus ( B 0 subscript 𝐵 0 B_{0} in GPa) for various ZB III-V semiconductors computed with OF-DFT/BLPS and the MGP and HC functionals with optimal choices of their parameters. The benchmark KS-DFT/BLPS results are also given. KS-DFT/BLPS and OF-DFT/BLPS with HC functional values are taken from Ref. 10 .

Systems	Functional	Parameters	$V_0(a.u/Cell)$	$E_0\left(eV\right)$	$B_0\left(GPa\right)$
GaAs	KS		274.2	-235.799	75
	HC	$\lambda =0.0130,\beta =0.783$	275.3	-235.782	81
	MGP	$a=0.434,b=0.524$	275.2	-235.801	75
GaSb	KS		354.2	-209.697	56
	HC	$\lambda =0.0100,\beta =0.720$	355.5	-209.739	58
	MGP	$a=0.341,b=0.752$	360.5	-209.696	49
GaP	KS		254.0	-243.079	80
	HC	$\lambda =0.0100,\beta =0.791$	255.0	-243.057	87
	MGP	$a=0.405,b=0.394$	252.5	-243.077	82
AlAs	KS		294.3	-232.908	80
	HC	$\lambda =0.0125,\beta =0.825$	301.6	-232.912	76
	MGP	$a=0.396,b=0.424$	296.2	-232.903	76
AlP	KS		274.2	-240.182	90
	HC	$\lambda =0.0120,\beta =0.845$	272.8	-240.199	91
	MGP	$a=0.378,b=0.319$	273.0	-240.180	81
AlSb	KS		382.0	-206.606	60
	HC	$\lambda =0.0120,\beta =0.750$	377.9	-206.588	61
	MGP	$a=0.359,b=0.832$	383.9	-206.606	55
InP	KS		310.7	-235.722	73
	HC	$\lambda =0.0120,\beta =0.885$	309.4	-235.696	66
	MGP	$a=0.394,b=0.351$	308.2	-235.724	62
InAs	KS		331.5	-228.537	65
	HC	$\lambda =0.0142,\beta =0.875$	334.0	-228.523	63
	MGP	$a=0.408,b=0.444$	330.1	-228.532	61
InSb	KS		424.5	-202.387	50
	HC	$\lambda =0.0120,\beta =0.810$	424.1	-202.381	49
	MGP	$a=0.361,b=0.850$	426.9	-202.386	46

Table 3. Table 3: Total computational wall time, number of truncated Newton (TN) optimization steps, total number of calls for evaluation of potential, wall time per call in evaluation of the total potential, and wall time per call in evaluation of kinetic potential with the WGC, MGP, and WT KEDFs. The wall time for evaluation of MGP kernel are also listed.

	KEDFs	2 atoms	16 atoms	128 atoms	1024 atoms
Total Wall time [s]	WGC	8.7	87.1	830.4	6264.2
	MGP	8.9	76.4	648.4	5486.4
	WT	5.0	51.5	456.8	4342.3
Number of optimized steps	WGC	10	10	11	12
	MGP	8	8	10	11
	WT	7	8	9	11
Number of calls for evaluation potentials	WGC	248	253	280	282
	MGP	191	212	224	228
	WT	173	207	219	236
Wall time of total potential[s/call]	WGC	0.03	0.32	2.81	21.27
	MGP	0.01	0.16	1.42	12.64
	WT	0.02	0.23	1.97	17.3
Wall time of kinetic potential[s/call]	WGC	0.03	0.27	2.37	17.88
	MGP	0.01	0.12	1.02	9.00
	WT	0.02	0.19	1.55	13.71
Wall time of evaluation kernel[s]	MGP	4.9	37.8	299.9	2383.3

Equations

$$E[\rho ]=T_{\mathrm{s}}[\rho ]+E_{\mathrm{H}}[\rho ]+E_{\mathrm{xc}}[\rho ]+E_{\mathrm{eN}}[\rho ]+E_{\mathrm{NN}},$$

$$T_{\mathrm{s}}[\rho ]$$

$$E_{\mathrm{H}}[\rho ]$$

$$E_{\mathrm{xc}}[\rho ]$$

$$\mathcal{L}[\rho ]=E[\rho ]-\mu \left(\int \rho (\mathbf{r})\mathrm{d}\mathbf{r}-N_e\right),$$

$$-\frac{1}{2}{\nabla }^2{\varphi }_i[\rho ](\mathbf{r})+\left[\underset{v_{\mathrm{s}}(\mathbf{r})}{\underbrace{v_{\mathrm{H}}[\rho ](\mathbf{r})+v_{\mathrm{xc}}[\rho ](\mathbf{r})+v_{\mathrm{ext}}(\mathbf{r})}}\right]{\varphi }_i[\rho ](\mathbf{r})={\epsilon }_i{\varphi }_i[\rho ](\mathbf{r}).$$

$$T\left[\rho \right]=T_{\mathrm{TF}}\left[\rho \right]+T_{\mathrm{vW}}\left[\rho \right]+T_{\mathrm{NL}}\left[\rho \right],$$

$$T_{\mathrm{NL}}\left[\rho \right]=\int {\rho }^{\alpha }\left(\mathbf{r}\right){\omega }_{T_{NL}}\left[\rho \right]\left(\mathbf{r},{\mathbf{r}}^{\prime }\right){\rho }^{\beta }\left({\mathbf{r}}^{\prime }\right)\mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime },$$

$$v_{{\mathrm{T}}_{\mathrm{s}}}[\rho ](\mathbf{r})+v_{\mathrm{H}}[\rho ](\mathbf{r})+v_{\mathrm{xc}}[\rho ](\mathbf{r})+v_{\mathrm{ext}}(\mathbf{r})=\mu .$$

$$4\pi {\rho }_{{\mathrm{T}}_{\mathrm{s}}}(\mathbf{r})=-4\pi \rho (\mathbf{r})-4\pi {\rho }_{\mathrm{xc}}(\mathbf{r})-4\pi {\rho }_N(\mathbf{r}),$$

$$\int {\rho }_{{\mathrm{T}}_{\mathrm{s}}}(\mathbf{r})\mathrm{d}\mathbf{r}=-N_e+1+N_N=1+\text{total charge},$$

$$E[\rho ]={\int }_0^{1}dt\int v[{\rho }_t](\mathbf{r})\frac{\mathrm{d}{\rho }_t(\mathbf{r})}{\mathrm{d}t}\mathrm{d}\mathbf{r},$$

$${\rho }_t(\mathbf{r})=t\rho (\mathbf{r})\to \frac{\mathrm{d}{\rho }_t}{\mathrm{d}t}=\rho (\mathbf{r}).$$

$$E[\rho ]=\int \mathrm{d}\mathbf{r}{\int }_0^{1}dtv[{\rho }_t](\mathbf{r})\rho (\mathbf{r}).$$

$$ϵ[\rho ](r)={\int }_0^{1}dtv[{\rho }_t](\mathbf{r}).$$

$$v[\rho ](\mathbf{r})=\int \mathrm{d}{\mathbf{r}}^{\prime }{\int }_0^1\mathrm{d}t\frac{\delta v[\rho ](\mathbf{r})}{\delta {\rho }_t({\mathbf{r}}^{\prime })}\rho ({\mathbf{r}}^{\prime })=\int \mathrm{d}{\mathbf{r}}^{\prime }{\int }_0^1\mathrm{d}t\frac{{\delta }^{2}E}{\delta \rho (\mathbf{r})\delta {\rho }_t({\mathbf{r}}^{\prime })}\rho ({\mathbf{r}}^{\prime }).$$

$$E[\rho ]=\int \mathrm{d}\mathbf{r}\int \mathrm{d}{\mathbf{r}}^{\prime }{\int }_0^1\mathrm{d}t{\int }_0^1\mathrm{d}{t}^{\prime }\frac{{\delta }^{2}E}{\delta {\rho }_t(\mathbf{r})\delta {\rho }_{t^{\prime }}({\mathbf{r}}^{\prime })}\rho (\mathbf{r})\rho ({\mathbf{r}}^{\prime }).$$

$$\frac{{\delta }^2{T}_{\mathrm{s}}}{\delta \rho (\mathbf{r})\delta \rho ({\mathbf{r}}^{\prime })}=-\frac{\delta v_s(\mathbf{r})}{\delta \rho ({\mathbf{r}}^{\prime })},$$

$${\chi }_{\mathrm{Lind}}(\eta )=-\frac{k_F}{{\pi }^2}\left(\frac{1}{2}+\frac{1-{\eta }^2}{4\eta }\ln \left|\frac{1+\eta }{1-\eta }\right|\right),$$

$$\hat{F}\bigg{[}\frac{\delta^{2}T_{\rm s}}{\delta\rho({\mathbf{r}})\delta\rho({\mathbf{r}^{\prime}})}\bigg{]}_{\rho({\mathbf{r}})=\rho_{0}}=-\frac{1}{\chi_{\text{Lind}}(\eta)}=\frac{\pi^{2}}{k_{F}}G_{\mathrm{Lind}}(\eta),$$

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

$$\frac{\delta T_{NL}[\rho ]}{\delta \rho (\mathbf{r})}=v_{{\mathrm{T}}_{\mathrm{NL}}}[\rho ](\mathbf{r})=\int \mathrm{d}{\mathbf{r}}^{\prime }{\int }_0^1\mathrm{d}t\frac{{\delta }^2{T}_{NL}}{\delta \rho (\mathbf{r})\delta {\rho }_t({\mathbf{r}}^{\prime })}\rho ({\mathbf{r}}^{\prime }),$$

$$\frac{{\delta }^2{T}_{NL}}{\delta \rho (\mathbf{r})\delta {\rho }_t({\mathbf{r}}^{\prime })}\simeq \frac{{\pi }^2}{(3{\pi }^2{)}^{\frac{1}{3}}}{\rho }^{-\frac{1}{6}}(\mathbf{r})G_{NL}[t{\rho }_0](|\mathbf{r}-{\mathbf{r}}^{\prime }|){\left(t\rho \right)}^{-\frac{1}{6}}({\mathbf{r}}^{\prime }),$$

$$G_{NL}(\eta )=G_{\mathrm{Lind}}(\eta )-3{\eta }^2-1.$$

$$v_{T_{NL}}[\rho ](\mathbf{r})=\frac{{\pi }^2}{(3{\pi }^2{)}^{1/3}}{\rho }^{-\frac{1}{6}}(\mathbf{r})\int d{\mathbf{r}}^{\prime }{\rho }^{\frac{5}{6}}({\mathbf{r}}^{\prime }){\int }_0^1\mathrm{d}t{t}^{-\frac{1}{6}}G[t{\rho }_0](\mathbf{r},{\mathbf{r}}^{\prime }).$$

$$v_{{\mathrm{T}}_{\mathrm{NL}}}(\mathbf{r})={\rho }^{-1/6}(\mathbf{r}){\hat{F}}^{-1}\left[\hat{F}\left[{\rho }^{5/6}(\mathbf{r})\right](q)\cdot {\omega }_T(q)\right].$$

$${\omega }_T(q)$$

$$={\omega }_{WT}(q)-\frac{6}{5}\frac{{\pi }^2}{(3{\pi }^2{)}^{1/3}}{\int }_0^{1}dt{t}^{\frac{5}{6}}\frac{d{G}_{NL}(\eta (q,t))}{dt}.$$

$${\tilde{v}}_{T_s}(q)=\frac{4\pi }{q^2}{\tilde{\rho }}_{T_s}(q)\to {\tilde{\rho }}_{T_s}(q)=\frac{1}{4\pi }{q}^2{\tilde{v}}_{T_s}(q),$$

$$\underset{q\to 0}{\lim }{\tilde{v}}_{T_s}(q)=\text{Constant},$$

$${\omega }_{MGP}(q)={\omega }_{T_{NL}}(q)+{\mathrm{erf}}^2(q)\frac{4\pi a}{q^2}{e}^{-b{q}^2}.$$

$$v_{T_s}(\mathbf{r})=$$

Full text

I Introduction

I.1 Orbital-free and Kohn–Sham DFT

Since the seminal work of Kohn and Sham Kohn and Sham (1965), the energy functional of the electron density can be written in terms of the noninteracting Kinetic energy (KE) functional, $T_{\rm s}[\rho]$, and the associated exchange-correlation (xc) functional, $E_{\rm xc}[\rho]$, as follows

[TABLE]

The nuclear-nuclear interaction, $E_{\rm NN}$ is density independent. The electron-nuclear interaction $E_{\rm eN}[\rho]=\int v_{\text{ext}}({\mathbf{r}})\rho({\mathbf{r}}){\mathrm{d}}{\mathbf{r}}$, is linear in the density, and the other terms [for closed-shell systems and introducing the Kohn–Sham orbitals, $\{\phi_{i}[\rho]({\mathbf{r}})\}$, including the yet unknown exchange–correlation (xc) energy density, $\varepsilon_{\rm xc}[\rho]({\mathbf{r}})$], take the form

[TABLE]

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

[TABLE]

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

[TABLE]

I.2 OF-DFT: The Kinetic energy functional conundrum

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

Local and semilocal (i.e., dependent on the value of the electron density and its gradient at a point in space) parametrizations of KEDFs have shown potential Karasiev and Trickey (2015); Trickey, Karasiev, and Jones (2009); Karasiev et al. (2009) and display a more than satisfactory agreement with KS-DFT for simulations of warm dense matter Karasiev et al. (2013). Recent work by the Della Sala group Laricchia et al. (2014) advocates for employing the Laplacian in addition to the gradient in the formulation of $T_{\rm s}[\rho]$ , and the possibility of reproducing accurate Kinetic energy densities Smiga et al. (2017).

Nonlocal versions of $T_{\rm s}[\rho]$, such as the Wang-Govind-Carter (WGC) Wang, Govind, and Carter (1998), the one proposed by Perrot Perrot (1994), the latest Huang-Carter (HC)Huang and Carter (2010), and others Chai and Weeks (2007); typically improve over local and semilocal functionals Xia and Carter (2012); Xia et al. (2012); Huang and Carter (2010); Ho, Lignères, and Carter (2008); Zhou, Lignères, and Carter (2005); Carling and Carter (2003) and all share the following form

[TABLE]

where $T_{\rm TF}\left[\rho\right]$ is the Thomas-Fermi Kinetic Energy Fermi (1927); Thomas (1927), $T_{\rm vW}\left[\rho\right]$ is the gradient correction due to von Weizsäcker von Weizsäcker (1935), and $T_{\rm NL}\left[\rho\right]$ is the remaining contribution which is by construction of nonlocal character. The general form of $T_{NL}$ is

[TABLE]

where the function $\omega_{T_{NL}}$ is commonly called “kernel of the nonlocal KEDF”, and $\alpha$ and $\beta$ are suitable exponents.

Unfortunately, the current state of the art for OF-DFT is that, although the algorithms have dramatically improved over the years and are typically orders of magnitude faster than KS-DFT, their application to semiconductors and molecular systems is still problematic Zhou, Lignères, and Carter (2005); Huang and Carter (2010); Xia and Carter (2012); Xia et al. (2012). For example, WGC Wang, Govind, and Carter (1998, 1999); Zhou, Lignères, and Carter (2005) and WT Wang and Teter (1992) KEDFs can reproduce remarkably accurate bulk properties for light metallic systems, but are less accurate for semiconductors and nonsimple metals. The HC KEDF Huang and Carter (2010) produces generally reliable bulk properties for semiconductors. However, it underestimates the electron density in bonding regions of CD Silicon and III-V semiconductors, it is less accurate than WT for metals, and it is substantially more computationally expensive than any other KEDFs.Huang and Carter (2010); Shin and Carter (2014) Density decomposition schemes that combine nonlocal with other functionals have found applicability to molecules, materials, and alloys Xia and Carter (2012); Shin and Carter (2014); Xia and Carter (2016) and thus appear to be an interesting avenue of research.

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

II Theoretical background

It is known that the exact $T_{s}$ functional satisfies a number of mathematical relations also known as exact conditions. It is arguable Perdew et al. (2005) that the more exact conditions a functional satisfies, the closer that functional necessarily is to the exact functional. It is through imposition of exact conditions that the Trickey group developed GGA KE functionals approaching the accuracy of KS-DFT in the limit of high temperature Karasiev et al. (2013). By imposing exact conditions on the asymptotic behavior of the linear response function, the Carter group formulated the HC KEDF which could be applied to systems with record inhomogenous density for OF-DFT simulations Huang and Carter (2010).

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

II.1 The Kinetic electron

The first unconstrained functional derivative of the DFT Lagrangian with respect to the density yields the Euler equation of DFT Liu and Ayers (2004). Namely,

[TABLE]

By taking the negative of the Laplacian of each term and rearranging, $-\nabla^{2}v_{\rm T_{s}}[\rho]({\mathbf{r}})=\nabla^{2}v_{\rm H}[\rho]({\mathbf{r}})+\nabla^{2}v_{\rm xc}[\rho]({\mathbf{r}})+\nabla^{2}v_{\rm ext}({\mathbf{r}})$. After solving the associated Poisson equations for each term in the long-range limit van Leeuwen and Baerends (1994), we find,

[TABLE]

where $\rho_{N}$ is the charge density of the nuclei (ions) and $\rho_{\rm xc}$ is related to the exchange-correlation hole. At long ranges, $\rho_{\rm xc}$ is identical to the xc hole van Leeuwen and Baerends (1994). Thus, at long ranges, $\rho_{\rm T_{\rm s}}$ is the Kinetic energy equivalent of the xc hole. However, upon simplification of the $4\pi$ terms and integration over all space, we find that

[TABLE]

Where the integration of the negative of the xc hole is 1, of the electron density is $N_{e}$, and of the nuclear density is $N_{N}$, and the total charge of the system is defined as $N_{N}-N_{e}$.

We define $\rho_{\rm T_{\rm s}}$ as the Kinetic electron, as for neutral systems it integrates to $+1$, in contrast to the xc hole which integrates to $-1$.

II.2 Hypercorrelation via a line integral

It is possible to recover an entire functional from the corresponding functional derivative using the following line integral van Leeuwen and Baerends (1995); Gaiduk, Chulkov, and Staroverov (2009):

[TABLE]

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

As the potential is the functional derivative $v({\mathbf{r}})=\frac{\delta E[\rho]}{\delta\rho({\mathbf{r}})}$, the result of the integration is path-independent. The scaled density, $\rho_{t}$, can then have almost any form van Leeuwen and Baerends (1995), the simplest of which is given by a linear interpolation between the vacuum and the density,

[TABLE]

Thus, the $t$-integral equation reduces to,

[TABLE]

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

[TABLE]

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

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

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

[TABLE]

It is important to notice that in the definition of the scaled second functional derivative on the rhs, $\rho_{t}$ appears only on one of the two functional derivatives.

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

[TABLE]

We should remark that due to the derivative discontinuity Perdew et al. (1982), the line integration in the above equations will cross discontinuities in the integrands when integrating the kernel to yield a potential, and also when integrating the potential to obtain an energy value. These discontinuities are integrable and do not pose formal problems.

II.3 FEG response

An exact condition for a KE functional is its relationship to the linear response function of the FEG. There is an established relation between the KE functional and the linear response function Wang and Carter (2000); Pearson, Smargiassi, and Madden (1993); Hohenberg and Kohn (1964). Taking a functional derivative of each term in the Euler equation Eq. (9),

[TABLE]

and the response function is given by \chi_{s}({\mathbf{r}},{\mathbf{r}^{\prime}})=\bigg{[}\frac{\delta v_{s}({\mathbf{r}})}{\delta\rho({\mathbf{r}^{\prime}})}\bigg{]}^{-1}. Unfortunately, the above relationship is only valid at self-consistency and trying to impose it directly would lead to impractical algorithms. However, for the FEG we have a simplified relationship, $\chi_{s}({\mathbf{r}},{\mathbf{r}^{\prime}})=\chi_{s}(|{\mathbf{r}}-{\mathbf{r}^{\prime}}|)$, resulting in a one-dimensional Fourier transform, $\chi_{\text{Lind}}(\eta)$, the Lindhard function Jones and Young (1971); Lindhard (1954),

[TABLE]

with $\eta=\frac{q}{2k_{F}}$ and the Fermi wavevector is given by, $k_{F}=(3\pi^{2}\rho_{0})^{1/3}$, with $\rho_{0}=\frac{N_{e}}{V}$.

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

[TABLE]

where $\hat{F}\left[\cdot\right]$ stands for Fourier transform in $q$, the conjugate variable of $|{\mathbf{r}}-{\mathbf{r}^{\prime}}|$. This exact condition is at the core of the original Hohenberg and Kohn functional Hohenberg and Kohn (1964) and inspired the entire class of nonlocal KE functionals, nicely summarized in Ref. 31. The function $G_{\mathrm{Lind}}(\eta)$ is introduced here,

[TABLE]

III Kinetic energy functional by integration

III.1 Imposing the response of the Free Electron Gas through hypercorrelation

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

[TABLE]

where $\frac{\delta v_{\rm T_{NL}}[\rho_{t}]({\mathbf{r}})}{\delta\rho_{t}({\mathbf{r}^{\prime}})}=\frac{\delta^{2}T_{NL}[\rho_{t}]}{\delta\rho({\mathbf{r}})\delta\rho_{t}({\mathbf{r}^{\prime}})}$ will be dependent on $\rho_{t}$ and consequently on $t$. Applying a local density approximation (LDA) on the polynomial terms of Eq. (20), $\frac{\delta^{2}T_{NL}}{\delta\rho({\mathbf{r}})\delta\rho_{t}({\mathbf{r}^{\prime}})}$ can be rewritten as:

[TABLE]

where, following the prescription of Eq. (7), and Eq. (19–20), the Fourier transform of $G_{NL}(|{\mathbf{r}}-{\mathbf{r}^{\prime}}|)$ is taken to be

[TABLE]

In the above, we have substituted $\rho_{0}\rightarrow t\rho_{0}$, which gives the following Fermi wavevector $k_{F}(t)=t^{\frac{1}{3}}\left(3\pi^{2}\rho_{0}\right)^{\frac{1}{3}}$, and the reciprocal space variable is $\eta({q},t)=\frac{q}{2k_{F}(t)}=\frac{q}{2t^{\frac{1}{3}}\left(3\pi^{2}\rho_{0}\right)^{\frac{1}{3}}}$. Thus, substituting Eq. (23) into the Eq. (22), the nonlocal Kinetic potential,$v_{T_{NL}}[\rho]({\mathbf{r}})$, can be rewritten as:

[TABLE]

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

Finally, the nonlocal Kinetic potential can be obtained by

[TABLE]

Where the kernel in reciprocal $\omega_{T}(q)$ is obtained from integration by parts

[TABLE]

The new kernel, $\omega_{T_{NL}}$, is given by the kernel of the WT KEDF Wang and Teter (1992) plus a correction term. We should point out that this correction term derives from assuming that the inverse Lindhard function is dependent on the electron density (this allowed us include in the line integral the average density $\rho_{0}$). Thus, we expect $\omega_{T_{NL}}$ in Eq. (III.1) to be most appropriate for substantially nonuniform densities. For light (simple) metallic systems, the distribution of electron density is very close to the uniform electron gas, thus WT kernel is already very good for these systems. We will see in our pilot calculations that for more complex systems, such as (nonsimple) metallic phases of Silicon, IV and III-V semiconductors, the correction term in conjuction with an additional correction arising from the presence of the “Kinetic electron” remarkably improves upon the performance of WT.

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

III.2 Imposing a nonzero Kinetic electron

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

[TABLE]

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

[TABLE]

and with that the Kinetic electron is absent.

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

[TABLE]

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

To understand the effect of the Kinetic electron, let us consider a very narrow Gaussian ($b\to+\infty$). This results in a spatially extended real-space Kinetic electron. This is the preferred shape for materials with a finite gap in which the dielectric screening is small. The opposite case ($b\to 0$) is preferred for metals, as there the dielectric screening is large and the real-space spatial extension of the Kinetic electron is reduced. Thus, we expect to need smaller $b$ values to model metallic systems, and large $b$ values to model semiconductors and insulators. These observations are consistent with the finding Pick, Cohen, and Martin (1970); Adler (1962); Kim and Görling (2002); Nazarov and Vignale (2011) (also exploited by Huang and Carter Huang and Carter (2010)) that, for $q\to 0$, the static response function for systems with gap, $\chi\to-\text{Constant}\times{q^{2}}$. An in-depth analysis of the relation between the low-$q$ behavior of $\chi$ and the long-range properties of the potential will be the subject of a follow-up work.

III.3 Introducing the MGP functional: details of the implementation

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

[TABLE]

[TABLE]

where $\hat{F}^{-1}\left[\right]$ indicates the inverse Fourier transform, and $v_{TF}({\mathbf{r}})+v_{vW}({\mathbf{r}})$ is the local part of the Kinetic energy potential composed of the sum of Thomas–Fermi and von Weizsäker terms. To avoid spurious dependencies from the parameter $a$, we choose such boundary conditions as $\omega_{MGP}(q=0)=0$.

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

[TABLE]

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

We mention here that the so-called LQ and HQ functionals Chai and Weeks (2007) (their potentials are constructed ad hoc to exactly reproduce the high-$q$ or the low-$q$ limits) feature a nonlocal kinetic energy that is also evaluated with a line integral similar to Eq. (32). The line integral in Ref. 11, however, is carried out only for the nonhomogeneous component of the density (i.e., $\rho({\mathbf{r}})-\rho_{0}$).

IV Computational details

The OF-DFT calculations are carried out with both a modified version of ATLASMi et al. (2016a) and PROFESS 3.0Chen et al. (2015). Calculations involving the HC KEDF are performed with PROFESS 3.0, all others are carried out in ATLAS and double checked with PROFESS 3.0 for consistency. The benchmark KS-DFT calculations are carried out with the same local pseudopotentials as the OF-DFT simulations. We employed the optimal effective local pseudopotentials (OEPP) Mi et al. (2016b) and the bulk-derived local pseudopotentials (BLPS) of Huang & Carter Huang and Carter (2008). KS-DFT calculations with OEPP are performed with Quantum-ESPRESSO (QE) Giannozzi et al. (2009), while the ones employing BLPS are carried out with ABINIT Gonze et al. (2002). In all calculations, the LDA xc functional by Perdew and Zunger Perdew and Zunger (1981) is adopted.

Crystal diamond (CD), body-centered cubic (BCC) and face-centered cubic (FCC) Silicon phases, as well as nine III-V ZB semiconductors are selected as benchmark systems. We compute total energy curves as a function of the lattice constant and extract the cell volume, $V_{0}$, the minimum energy, $E_{0}$ and the bulk modulus, $B_{0}$. These were computed using the prescription of Carter and coworkers Huang and Carter (2010) fitting the energy curves vs volume against Murnaghan’s equation of state Murnaghan (1944) within 0.95 $V_{0}$ to 1.05 $V_{0}$.

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

V Results and Discussion

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

V.1 Bulk properties for both metallic and semiconductor phases of Silicon

The optimal parameters for MGP functionals along with the optimal $\gamma$ and $\rho^{*}$ for WGC, optimal $\beta$ and $\lambda$ for HC functionals and the calculated quantities for all KEDFs are collected in Table 1.

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

Modelling CD Si with OF-DFT has historically been a challenge which has been addressed by several nonlocal functionals with density-dependent kernel Wang, Govind, and Carter (1998, 1999); Shin and Carter (2014); Huang and Carter (2010). For example Zhou, Lignères, and Carter (2005), WT is unable to reproduce a bound curve for CD Si. As MGP functional has a density-independent kernel, we had no expectations that CD Si would be modeled correctly. Table 1, however, shows that MGP is capable of producing bound energy curves and overall bulk properties that are close to the KS-DFT benchmark compared to other KEDFs. To the best of our knowledge, MGP is the only KEDF with density-independent kernel capable of simultaneously reproducing KS-DFT equilibrium total energies, bulk moduli, and equilibrium volumes for both metallic and semiconductor phases of Silicon.

We also tested a modified version of MGP only including WT and Kinetic electron term. We found that the modified functional could not reproduce the benchmark despite efforts in optimizing the $a$ and $b$ parameters. This indicates that, to achieve accurate results it is not enough to simply add the correction terms separately. Instead, the correction term derived from functional integration and the Kinetic electron must be included together for MGP to approach both metallic and semiconducting phases.

V.2 Benchmarks for III-V semiconductors

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

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

V.3 Analysis of the electron density

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

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

V.4 Computational efficiency

To evaluate the computational efficiency of MGP, we compare the computational cost for optimizing the electron density in CD Si (2 atoms), $2\times 2\times 2$ (16 atoms), $4\times 4\times 4$ (128 atoms), and $8\times 8\times 8$ (1024 atoms) supercells with WGC, WT, and MGP KEDFs. All calculations are performed on single thread with PROFESS employing the BLPS pseudopotentials using a large Kinetic energy cutoff (4000 eV). For all KEDFs we employ optimal parameters (WGC: $\gamma=4.2$, $\rho^{*}$=$\rho_{0}$; WT: $\alpha=\beta=5/6$; MGP: $a$=0.364 and $b$=0.570 ). Since the computational cost of HC KEDF is much higher than that of other functionals (about 1000 times higher than MGP for CD Si) Huang and Carter (2010); Shin and Carter (2014), we exclude it in the comparison. In principle, the only difference in computational cost between MGP and WT KEDFs is the initial kernel building step in reciprocal space.

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

Among WGC, MGP and WT KEDFs, WGC is the most expensive KEDF in terms of total wall time. Additionally, we find that MGP KEDFs converges in all cases considered with 12 iterations of a truncated Newton minimization. In contrast, WGC KEDF is numerical unstable in some cases (this is a well-known limitation of this functional Huang and Carter (2010); Shin and Carter (2014)).

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

VI Conclusions

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

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

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

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

Our analysis shows that there is room for improvement, particularly for two aspects of the MGP functional. First, the Kinetic electron can be better parametrized following prescriptions that have been very successful in formulating exchange functionals in real space Perdew (1985); Gritsenko et al. (1995). Secondly, there exist formulations of “jellium with gap” models for the Lindhard function Constantin et al. (2017) which have been recently (and successfully) applied to the formulation of GGA KEDF Smiga et al. (2017). These can be extended to the line-integral formulation of MGP. Thirdly, a density-dependent kernel version of MGP can be constructed to satisfy additional exact conditions, such as density and coordinate scaling relations Levy (1988).

VII Supplementary Material

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

VIII Acknowledgments

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