Orbital-Free DFT Correctly Models Quantum Dots When Asymptotics, Nonlocality and Nonhomogeneity Are Accounted For

TL;DR

This paper develops a new family of nonlocal kinetic energy functionals for orbital-free DFT, enabling accurate simulations of quantum dots and clusters with highly inhomogeneous electron densities, previously a major limitation.

Contribution

The authors introduce a universal, nonlocal kinetic energy functional that correctly models asymptotics and nonhomogeneity, extending OF-DFT's applicability to quantum dots and clusters.

Findings

Achieves close to chemical accuracy in electronic energy calculations.

Reproduces electron density within 5% of benchmark results.

Enables orbital-free DFT to simulate highly inhomogeneous systems.

Abstract

Million-atom quantum simulations are in principle feasible with Orbital-Free Density Functional Theory (OF-DFT) because the algorithms only require simple functional minimizations with respect to the electron density function. In this context, OF-DFT has been useful for simulations of warm dense matter, plasma, cold metals and alloys. Unfortunately, systems as important as quantum dots and clusters (having highly inhomogeneous electron densities) still fall outside OF-DFT's range of applicability. In this work, we address this century old problem by devising and implementing an accurate, transferable and universal family of nonlocal Kinetic Energy density functionals that feature correct asymptotics and can handle highly inhomogenous electron densities. For the first time to date, we show that OF-DFT achieves close to chemical accuracy for the electronic energy and reproduces the electron density to about 5\% of the benchmark for semiconductor quantum dots and metal clusters. Therefore, this work demonstrates that OF-DFT is no longer limited to simulations of systems with nearly homogeneous electron density but it can venture into simulations of clusters and quantum dots with applicability to rational design of novel materials.

Tables

Table 1. Table 1: Mean unsigned error (MUE) for the total energy compared to KS-DFT in eV/atom. Superscripts S 𝑆 S and V ​ S 𝑉 𝑆 VS stand for “strained” and “very strained”, respectively.

Systems	LMGP	LMGP0	LWT	TF+$\frac{1}{5}$vW	WT
Mg8	0.18	0.63	1.19	1.09	8.79
Si8	0.22	2.17	4.86	1.46	41.7
Ga4As4	0.34	2.21	6.15	1.55	51.8
Mg50	0.05	0.35	0.84	0.95	3.23
Si50	0.11	0.95	3.73	1.53	16.4
Ga25As25	0.13	1.06	4.29	1.67	22.7
Mg${}_{}{}^S{}_8{}^{}$	0.28	1.16	2.66	0.27	19.4
Mg${}_{}{}^{VS}{}_8{}^{}$	0.09	1.67	3.88	0.10	24.0

Table 2. Table 2: Mean unsigned relative error (MURE) for the electron density (measured by this quantity, 1 2 ​ N e ​ ∫ | ρ OF-DFT ​ ( 𝐫 ) − ρ KS-DFT ​ ( 𝐫 ) | ​ 𝑑 𝐫 1 2 subscript 𝑁 𝑒 subscript 𝜌 OF-DFT 𝐫 subscript 𝜌 KS-DFT 𝐫 differential-d 𝐫 \frac{1}{2N_{e}}\int|\rho_{\text{OF-DFT}}({\mathbf{r}})-\rho_{\text{KS-DFT}}({\mathbf{r}})|d{\mathbf{r}} ) in percentage points. Superscripts S 𝑆 S and V ​ S 𝑉 𝑆 VS stand for “strained” and “very strained”, respectively.

Systems	LMGP	LMGP0	LWT	TF+$\frac{1}{5}$vW	WT
Mg8	3.79	4.12	4.05	11.36	16.0
Si8	4.84	4.90	4.74	8.28	17.5
Ga4As4	5.40	5.43	4.89	8.94	19.3
Mg50	3.31	3.42	2.38	9.56	10.3
Si50	4.59	4.65	3.60	7.24	14.2
Ga25As25	5.21	5.26	3.19	7.79	16.8
Mg${}_{}{}^S{}_8{}^{}$	5.20	5.34	5.29	7.63	18.6
Mg${}_{}{}^{VS}{}_8{}^{}$	3.94	4.10	4.87	5.60	17.5

Equations

$$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 }_{\mathrm{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{NL}}(\mathbf{r})=\rho (\mathbf{r}{)}^{-\frac{1}{6}}{\hat{F}}^{-1}\left[\hat{F}\left[{\rho }^{\frac{5}{6}}\right](q){\omega }_{MGP}(q)\right](\mathbf{r}),$$

$${\omega }_{MGP}(q)={\omega }_{WT}(q)+{\omega }_{Hyper}(q)+{\omega }_e(q).$$

$${\omega }_e(q)=\frac{4\pi a}{q^2}{\mathrm{erf}}^2(q)\exp (-a{q}^2),$$

$$\displaystyle v_{{\rm NL}}[\rho({\mathbf{r}})]({\mathbf{r}})\xleftarrow{~{}\rm spline~{}}\bigg{\{}v_{{\rm NL}}[\rho_{1}]({\mathbf{r}}),\ldots,v_{{\rm NL}}[\rho_{k}]({\mathbf{r}}),\ldots,v_{{\rm NL}}[\rho_{N_{sp}}]({\mathbf{r}})\bigg{\}}.$$

$$T_{\mathrm{NL}}[\rho ]=\int d\mathbf{r}\int dt\rho (\mathbf{r})v_{\mathrm{NL}}[t\rho ](\mathbf{r}).$$

Full text

Metal clusters and quantum dots constitute an important class of systems of pivotal importance for materials design particularly in photovoltaics Lan et al. (2014), catalysis Tyo and Vajda (2015), and even quantum computing Huo et al. (2013). Although these fields are already strongly shaped by computer-aided design, the high computational cost of available quantum-mechanical methods such as Kohn–Sham density functional theory (KS-DFT) Maitra (2016); Burke (2012) is hampering futher progress. In this playing field, what is really needed is a breakthrough in techniques alternative to the current standard, and among them Giese and York (2017); VandeVondele et al. (2012); Krishtal et al. (2015); Yang (1991); Jacob and Neugebauer (2014); Fedorov (2017); Goedecker (1999) Orbital-Free Density Functional Theory (OF-DFT) is a promising candidate.

OF-DFT is a promising and intriguing alternative because approximate density functionals for the kinetic energy entirely replace the need to solve for a Schrödinger equation. This completely bypasses its inherent complexity. Particularly, OF-DFT algorithms are promising because they involve a computational scaling of at most $\mathcal{O}(N\ln N)$, where $N$ is a measure of the system size, and a memory requirement of only $\mathcal{O}(N)$ Wesolowski and Wang (2013); Karasiev and Trickey (2012); Witt et al. (2018).

Unfortunately, even though OF-DFT has already proven to be successful for simulations of million-atom systems involving crystalline and liquid metals and alloys Witt et al. (2018); Shao et al. (2018); Hung and Carter (2009); González et al. (2004), as well as plasmas and warm dense matter White et al. (2013); Ding et al. (2018); Sjostrom and Daligault (2014), its applicability has been severely limited by the accuracy of the available Kinetic Energy density functionals (KEDFs). For example, finite systems such as metal clusters and quantum dots have been outside the range of applicability of OF-DFT.

In this work, we achieve a breakthrouth by carefully balancing three important aspects defining the KEDFs: asymptotics of the corresponding potential, intrinsic nonlocality, and ability to handle nonhomogeneous systems. Thus, already at conception, we make sure that the functionals are nonlocal, that their asymptotics matches the known exact behavior, and that their nonlocal kernels adapt to such large density inhomogeneities as the ones occurring at the interface of nonperiodic systems with the vacuum. For finite systems, such as clusters, the latter is perhaps the most important aspect, which we tackle head on.

In the following, we cast our KEDF development in the current stat-of-the-art and derive the main theoretical and algorithmmical developments. Afterwards, we benchmark the resulting density functionals by carrying out OF-DFT simulations on 4 metal clusters and 4 semiconductor quantum dots realized in one hundred possible geometries for each, spanning energy windows of up to several tens of eV.

Over the past two decades there have been tremendous advances in OF-DFT development. Various KEDFs have been proposed Constantin et al. (2018); Luo et al. (2018); Karasiev and Trickey (2015); Trickey et al. (2009); Karasiev et al. (2009, 2013); Wang et al. (1998, 1999); Perrot (1994); Laricchia et al. (2014); Smiga et al. (2017); Chai and Weeks (2007); Xia and Carter (2012); Xia et al. (2012); Huang and Carter (2010); Ho et al. (2008); Zhou et al. (2005); Carling and Carter (2003). The majority of these functionals are appropriate for main-group metallic bulk systems, and some show potential for modeling bulk semiconductors Huang and Carter (2010); Mi et al. (2018). Although semilocal KEDFs Karasiev and Trickey (2015); Trickey et al. (2009); Karasiev et al. (2009) have seen a recent resurgence Luo et al. (2018); Constantin et al. (2018), nonlocal KEDFs (such as WGC Wang et al. (1998, 1999), HC Huang and Carter (2010), WT Wang and Teter (1992), MGP Mi et al. (2018), and others Pearson et al. (1993); Smargiassi and Madden (1994); Perrot (1994)) have historically delivered better results, particularly for bulk solids. An inspiring study by Chan, Cohen and Handy found semilocal KEDFs to be theoretically appropriate for applications to clusters Chan et al. (2001). This was followed by several works on metallic clusters Aguado (2001); Aguado et al. (2001, 1999) which employed simple combinations of Thomas–Fermi (TF) and fractions of von Weizsäcker (vW) KEDFs (e.g., $\frac{1}{9}$ or $\frac{1}{5}$). Unfortunately, the conclusions of those studies were mixed. Thus, in this work we depart from semilocal KEDFs and adopt fully nonlocal ones exploiting the typical nonlocal KEDF ansatz,

[TABLE]

where $T_{\rm TF}$ is the Thomas–Fermi kinetic energy Fermi (1927); Thomas (1927), $T_{\rm vW}\left[\rho\right]$ is the von Weizsäcker functional von Weizsäcker (1935), and $T_{\rm NL}\left[\rho\right]$ is the remaining nonlocal contribution. The general form of $T_{\rm NL}$ is

[TABLE]

Where, $\omega_{\rm NL}\left[\rho\right]\left({\mathbf{r}},{\mathbf{r}^{\prime}}\right)$ is the kernel of the nonlocal KEDF, and $\alpha$ and $\beta$ are positive scalars.

Let us first describe details of the MGP functional and then outline the new developments that extend its applicability to finite systems. Following the lead of Kohn and Sham Kohn and Sham (1965), the starting point of our derivations is the linear response of the Free Electron Gas. This starting point is common among nonlocal functionals Wang and Carter (2000). After a procedure of functional integration, this yields the WT functional Wang and Teter (1992) in the zero-th order and the MGP functional at higher orders Mi et al. (2018). MGP’s kernel correctly behaves in the low $q$ limit by construction, as we impose the so-called “kinetic electron” (vide infra) and therefore it can potentially approach systems beyond bulk metals. In our previous work Mi et al. (2018), we found that MGP reproduces with remarkable accuracy the energetic properties and electron densities of silicon and several III–V semiconductors provided that two free parameters are adjusted.

MGP’s nonlocal potential is given by

[TABLE]

where the reciprocal space variable, $q$, corresponds to $|{\mathbf{r}}-{\mathbf{r}^{\prime}}|$, and $\hat{F}[\cdot]$ stands for Fourier transform. Thus, the inherent approximation in Eq. (3) is that the kernel only depends on the magnitude of $|{\mathbf{r}}-{\mathbf{r}^{\prime}}|$. The kernel introduced in Eq. (3), takes the following form

[TABLE]

The first term, $\omega_{WT}$, is the kernel of the WT functional. The second term, $\omega_{Hyper}$, originates from functional integration which, for historical reasons Hessler et al. (1999), we also call hypercorrelation Mi et al. (2018). The functional integration step is employed to transform the kernel borrowed from the response function of the free electron gas into a kernel that can be used in the computation of the KEDF potential. In practice, the integration is carried out numerically after an integration by parts (see additional details in Ref. 40). The third term, $\omega_{e}$, is the contribution encoding the kinetic electron and needs to be approximated. The kinetic electron arises from the long-range behavior of the exchange potential Mi et al. (2018). An exact condition for the exchange potential in finite systems is that the negative of its source (i.e., $\frac{1}{4\pi}\nabla^{2}v_{xc}(r)$) integrates to unity. As exchange potential and KEDF potential are connected via the Euler equation of DFT, the source of the KEDF (i.e., the kinetic electron) shall integrate to the opposite of the exchange hole.

In contrast to Ref. 40, here we propose a universal form for $\omega_{e}$, containing no adjustable parameters. Namely,

[TABLE]

where $a$ is a parameter that we relate to the number of electrons, $a=A/N_{e}^{2/3}$, with $A=0.2$ and $N_{e}$ is the total number of electrons.

Thus, if the kernel includes $\omega_{WT}(q)$, $\omega_{WT}(q)+\omega_{Hyper}(q)$ or $\omega_{WT}(q)+\omega_{Hyper}(q)+\omega_{e}(q)$ in Eq.(4), then the corresponding KEDF is called WT, MGP0, or MGP. These kernels are only dependent on the average electron density, $\rho_{0}$, through $\eta=q/(2k_{F})$, where $k_{F}=(3\pi^{2}\rho_{0})^{1/3}$ is the Fermi wave vector associated with the average density. This approximation is too strong and is the source reason for needing adjustable parameters in these functionals. Thus, the kernel should be made dependent on the total electron density, $\rho({\mathbf{r}})$, instead of $\rho_{0}$. In principle, this would benefit and improve the description of systems where the distribution of electron density strongly deviates from uniformity.

Such a strong approximation is shared among most nonlocal functionals. In this work we propose a method that tackles this limitation and in Figure 1 we hint at the proposed workaround.

Inspired by the success of the local density approximation (LDA) Ceperley and Alder (1980); Perdew and Wang (1992), we assume that the potential at a point ${\mathbf{r}}$ in space can be approximated by the one of a nonlocal functional evaluated with a kernel $\omega[\rho_{0}=\rho({\mathbf{r}})](q)$. This is the same principle as LDA applied to nonlocal kernels rather than to the energy densities, as customarily done. Unfortunately, implementing this principle directly would imply a super-quadratic computational cost ($N$ kernel evaluations each costing $N\ln N$). Fortunately, the computational scaling can be brought back to \mathcal{O}\big{(}N\ln N\big{)} by employing spline techniques. Figure 1 provides a visual depiction of the proposed spline method which we explain in detail in the following.

A series of $\rho_{0}$ values is considered, $\{\rho_{k}\}$, obtained by dividing the range between 0 and $\rho_{max}=\text{max}[\rho({\mathbf{r}})]$, in equally spaced segments and choosing the total number of $k$ points to be $N_{sp}=40$. For each $\rho_{k}$, there is a corresponding kernel, $\omega_{\rm NL}(q,\eta(\rho_{k}))$, that can be tabulated and recovered ahead of computing the potential. Thus, a series of nonlocal potentials is obtained based on Eq. (3), \big{\{}v_{{\rm NL}}[\rho_{k}]({\mathbf{r}})\big{\}}_{k\in 1\ldots N_{sp}}, and the LDA procedure can be exploited invoking splines,

[TABLE]

This is a scheme for constructing LDA versions of MGP (LMGP), WT (LWT) and MGP0 (LMGP0) functionals from the corresponding kernels.

Finally, The nonlocal contribution to the total kinetic energy is recovered by a second functional integration

[TABLE]

Two other prescriptions for generating density-dependent kernels exist. WGC Wang et al. (1999) exploits a Taylor expansion for the kernel around a reasonable value near $\rho_{0}$. Unfortunately, this can result in numerical instabilities when the electron density distribution differs much from the one of the free electron gas. Another example is the kernel of the HC functional Huang and Carter (2010). It is obtained by solving a differential equation when the electron density is updated. To ameliorate the computational cost, Huang and Carter offer an implementation of HC also featuring a spline technique in the spirit of Soler and coworkers Román-Pérez and Soler (2009). Despite this, the computational cost of HC compared to WGC is still orders of magnitude larger. In terms of performance, the WGC functional can reproduce well KS-DFT results for main group bulk metals; HC functional can achieve significant improvement over previous functionals for semiconductors, showing promise for simulating finite systems when its two free parameters ($\lambda$ and $\beta$) are appropriately adjusted Xia and Carter (2012).

A major advantage of the LWT and LMGP family of functionals compared to WGC and HC is the fact that they are universal functionals with no adjustable parameters. One issue with universal functionals is that they may not be transferable. I.e., they may work well for a certain system, but less well for others. In the following, we craft strict and conservative benchmarks for the proposed functionals that undeniably show their superiority compared to the current state-of-the-art and their transferability among an array of cluster sizes and types.

To assess the accuracy of the new family of KEDFs, we choose random clusters generated by CALYPSO Wang et al. (2012, 2010); Lv et al. (2012). Standard KS-DFT calculations provide benchmark values for the total energy and electron density (KS-DFT employs the exact KEDF) are carried out with Quantum ESPRESSO (QE) Giannozzi et al. (2009). The OF-DFT simulations are performed with a modified version of ATLAS Mi et al. (2016a); Shao et al. (2016) and PROFESS 3.0 Chen et al. (2015). To avoid dealing with nonlocal pseudopotentials, the optimal effective local pseudopotentials (OEPP) Mi et al. (2016b) are employed for both OF-DFT and KS-DFT calculations. LDA exchange-correlation energy functional by Perdew and Zunger Perdew and Zunger (1981) is employed in all calculations. Additional technical and computational details are available in the Supplementary Materials.

With KS-DFT quantities as reference, we initially select two types of clusters: Mg8 and Si8. For each, we compute the total energy of 100 random structures and collect the computed values in Figure 2. The figure shows a progressive improvement when adopting the functionals WT $\to$ MGP $\to$ LMGP. In particular, we note that the consistent bias of MGP (slope differing from KS-DFT) is eliminated by the LDA procedure in LMGP. As shown in the lower insets of Figure 2, MGP improves dramatically total energies in comparison with WT. Furthermore, the three new parameterless functionals (LWT, LMGP0, and LMGP) are found to outperform their $\rho_{0}$-dependent kernel counterparts. We should note that while LWT and LMGP/LMGP0 functionals are universal (no adjustable parameters), MGP results are obtained by adjusting the parameter $a$ independently for Mg8 and Si8 clusters ($a=0.35$ and $0.5$, respectively). Strikingly, LMGP values are found to be essentially on-top of the KS-DFT results, providing us with an indication that the LDA procedure implemented by splines is stable and accurate for these systems.

To confirm the transferability of our new functionals, we select six additional cluster systems: Mg50, Si50, Ga4As4, Ga25As25, and two more Mg8 (i.e., Mg${}_{8}^{S}$ and Mg${}_{8}^{VS}$) featuring shorter average interatomic distances. The new set of systems provides us with an array of metallic to semiconducting quantum dot-like clusters. As shown in Figure 3, the performance of our new functionals is consistent for all systems reproducing total energies across a wide window of energy spanning several eV per atom. In terms of absolute values of total energies, LWT and LMGP0 results are higher and lower than KS-DFT results, respectively. These results are source of considerable excitement – not only the LMGP energy values correlate almost perfectly with the KS-DFT benchmark, but more importantly the LDA procedure (which is numerical in nature) is found to be stable for all systems considered. LMGP converges for more than 90% of the structures in all systems with an average convergence rate of over 95%.

To quantify the performance of our new functionals, Table 1 shows the mean unsigned error (MUE) of total OF-DFT energies compared to KS-DFT for the 100 random structures of each system computed with LWT, LMGP, and LMGP0, as well as WT, and TF+$\frac{1}{5}$vW.

An option is to also compare against the Huang-Carter Huang and Carter (2010) (HC) and the Wang-Govind-Carter (WGC) functionals. However, even though HC is considered to be the most accurate KEDF presently available, it is also known for drawbacks that make it unsuitable for applications to finite systems Xia and Carter (2012). Xia and Carter Xia and Carter (2012) found that especially for molecules and solid surfaces, the computational cost of HC is hundreds of times higher than the WGC functional. Additionally, despite our best efforts, HC (with $\lambda>0$, an appropriate value for systems with gap, such as clusters) as well as WGC (second-order Taylor expansion of the kernel) did not converge for more than 10 of the 100 cluster structures considered for all systems. Thus, in this work, we compare against the other functionals.

Table 1 shows that LWT considerably improves over WT. This indicates that the LDA procedure improves significantly the corresponding nonlocal KEDF with density independent kernel while at the same time removing the $\rho_{0}$ dependence in the functional. Interestingly, LMGP0 performs even better than LWT. This is an indication that the hypercorrelation term in the kernel further improves the performance of the functional. Adding the kinetic electron (i.e., the additional term $\omega_{e}(q)$ in the kernel, see Eq. (4)), the LMGP functional achieves additional and important improvement over LMGP0, lining itself up to the KS-DFT results in an almost quantitative fashion. Strikingly, this is so despite the relatively uncomplicated formalism for the kinetic electron term in Eq. (5).

Reproducing the electron density is as important as reproducing the total energy. This was pointed out recently for exchange-correlation functionals Medvedev et al. (2017) and it is even more important for KEDFs. Thus, we define $\frac{1}{2N_{e}}\int|\rho_{\text{OF-DFT}}({\mathbf{r}})-\rho_{\text{KS-DFT}}({\mathbf{r}})|d{\mathbf{r}}$, a measure of the electron density difference between KS-DFT and OF-DFT. The performance of the various functionals in reproducing the KS-DFT electron density are listed in Table 2. Once again, the three new functionals constitute an improvement over TF+$\frac{1}{5}$vW and WT functionals. We point out that although the total energies computed by TF+$\frac{1}{5}$vW only partially differ from LWT, the LWT electron density is of much higher quality than the one from TF+$\frac{1}{5}$vW.

In conclusion, we have addressed a long standing problem in the field of OF-DFT, i.e., the simulation of finite systems, such as quantum dots, reproducing benchmark KS-DFT results with unprecedented accuracy. This constitutes a major step forward for OF-DFT, a framework that was thought to only be applicable to bulk metals and alloys. Quantum dot structure prediction is now feasible with OF-DFT, opening the door to a new regime of applicability for this method.

Our results are achieved by (1) imposing correct asymptotic behavior of the kinetic energy potentials, (2) accounting for nonlocality in the functional by construction, and (3) allowing the functional to adapt to systems with highly inhomogeneous electron density, via a technique inspired by the LDA approximation. Such a comprehensive, yet simple prescription leads to a numerically stable family of noninteracting kinetic energy functionals which we apply to 8 different quantum dots each realized in 100 different structures spanning energy windows ranging between 10 and 80 eV. Our most refined functional, LMGP, consistently reproduces the KS-DFT electronic energy for all 50-atom quantum dots to within 130 meV/atom. The energies of the 8-atom clusters are within 340 meV/atom of the KS-DFT reference. These errors are for the most part systematic in nature, as the OF-DFT energy values correlate almost perfectly to the KS-DFT benchmarks.

Although the noninteracting kinetic energy functionals presented here predict with unprecedented accuracy the total energy and electron density of the considered quantum dots, there is still room for improvement both in terms of computational accuracy as well as efficiency. Particularly, LMGP shows a significant improvement in comparison to LWT in terms of total energies, but struggles to improve the quality of the electron density. This indicates that the simple LDA approximation for the kernels can be further improved, for instance by including a dependency over the density gradient. This is currently being investigated.

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1553993.