Exact kinetic energy functional for general many-electron systems

TL;DR

This paper derives the exact kinetic energy functional for general many-electron systems using a two-density model, addressing a longstanding challenge in orbital-free density functional theory.

Contribution

It introduces a novel two-density model approach that yields the exact kinetic energy functional, improving the theoretical foundation of orbital-free DFT.

Findings

Exact kinetic energy functional derived for many-electron systems.

Functional includes previously suggested functionals, excluding Thomas-Fermi.

Advances orbital-free DFT by providing a precise functional form.

Abstract

The exact form of the kinetic energy functional has remained elusive in orbital-free models of density functional theory (DFT). This has been the main stumbling block for the development of a general-purpose framework on this basis. Here, we show that on the basis of a two-density model, which represents many-electron systems by mass density and spin density components, we can derive the exact form of such a functional. The exact functional is shown to contain previously suggested functionals to some extent, with the notable exception of the Thomas-Fermi kinetic energy functional.

Equations

$$n(\mathbf{r})=\rho (\mathbf{r})+S(\mathbf{r}).$$

$$\Psi (\mathbf{r})={\rho }^{\text{nicefrac}12}(\mathbf{r})+\mathbf{i}{\mathbf{e}}_S(\mathbf{r})S^{\text{nicefrac}12}(\mathbf{r}).$$

$$H(\mathbf{r})=-\frac{1}{2}{\nabla }^2+v_0(\mathbf{r})+\mathbf{i}{\mathbf{e}}_b(\mathbf{r})v_b(\mathbf{r})=:{\hat{h}}_0(\mathbf{r})+\mathbf{i}{\mathbf{e}}_b(\mathbf{r})v_b(\mathbf{r}).$$

$$\left(\begin{array}{cc}{\hat{h}}_0& -v_b\\ v_b& {\hat{h}}_0\end{array}\right)\left(\begin{array}{c}{\rho }^{\text{nicefrac}12}\\ S^{\text{nicefrac}12}\end{array}\right)=\mu \left(\begin{array}{c}{\rho }^{\text{nicefrac}12}\\ S^{\text{nicefrac}12}\end{array}\right),$$

$$v_b=\frac{1}{\rho +S}\left[S^{\text{nicefrac}12}\left(-\frac{{\nabla }^2}{2}\right){\rho }^{\text{nicefrac}12}-{\rho }^{\text{nicefrac}12}\left(-\frac{{\nabla }^2}{2}\right)S^{\text{nicefrac}12}\right].$$

$$\mathbf{i}{\mathbf{e}}_S(\mathbf{r})S^{\text{nicefrac}12}(\mathbf{r})=\sum _{j=1}^3\mathbf{i}{\mathbf{e}}_j{S}_j^{\text{nicefrac}12}(\mathbf{r}).$$

$$\mathbf{i}{\mathbf{e}}_b(\mathbf{r})v_b(\mathbf{r})=\sum _{j=1}^3\mathbf{i}{\mathbf{e}}_j{v}_j(\mathbf{r}).$$

$$\left(\begin{array}{cccc}{\hat{h}}_0& -v_1& -v_2& -v_3\\ v_1& {\hat{h}}_0& v_3& -v_2\\ v_2& -v_3& {\hat{h}}_0& v_1\\ v_3& v_2& -v_1& {\hat{h}}_0\end{array}\right)\left(\begin{array}{c}{\rho }^{\text{nicefrac}12}\\ S_1^{\text{nicefrac}12}\\ S_2^{\text{nicefrac}12}\\ S_3^{\text{nicefrac}12}\end{array}\right)=\mu \left(\begin{array}{c}{\rho }^{\text{nicefrac}12}\\ S_1^{\text{nicefrac}12}\\ S_2^{\text{nicefrac}12}\\ S_3^{\text{nicefrac}12}\end{array}\right).$$

$$v_k=\frac{1}{4}\frac{{\hat{\mathbf{J}}}_k-B_k^j{\hat{\mathbf{J}}}_i-B_k^i{\hat{\mathbf{J}}}_j}{1-B_k^j{\chi }_i^{+}-B_k^i{\chi }_j^{-}},$$

$${\hat{\mathbf{J}}}_k$$

$$\nabla \left({\rho }^{\text{nicefrac}12}\nabla S_k^{\text{nicefrac}12}-S_k^{\text{nicefrac}12}\nabla {\rho }^{\text{nicefrac}12}+S_i^{\text{nicefrac}12}\nabla S_j^{\text{nicefrac}12}-S_j^{\text{nicefrac}12}\nabla S_j^{\text{nicefrac}12}\right)$$

$$B_k^j$$

$${\chi }_k^{\pm }$$

$${ϵ}_{\text{tot}}={\rho }^{\text{nicefrac}12}{\hat{h}}_0{\rho }^{\text{nicefrac}12}+\sum _k{S}_k^{\text{nicefrac}12}{\hat{h}}_0{S}_k^{\text{nicefrac}12},$$

$${\epsilon }_{\text{e-n}}\left(\mathbf{r}\right)=v_{\text{e-n}}\left(\mathbf{r}\right)n\left(\mathbf{r}\right),$$

$${\epsilon }_{\text{e-e}}\left(\mathbf{r}\right)=v_{\text{e-e}}\left(\mathbf{r}\right)n\left(\mathbf{r}\right),$$

$${\epsilon }_{\text{k}}(\mathbf{r})=-\frac{1}{4}{\nabla }^{2}n(\mathbf{r})+\frac{1}{8}\left(\frac{{\left[\nabla \rho (\mathbf{r})\right]}^2}{\rho (\mathbf{r})}+\sum _{k=1}^3\frac{{\left[\nabla S_k(\mathbf{r})\right]}^2}{S_k(\mathbf{r})}\right).$$

$${\epsilon }_{\text{k}}(\mathbf{r})={\epsilon }_{\text{k}}^1(\mathbf{r})+{\epsilon }_{\text{k}}^2(\mathbf{r})+{\epsilon }_{\text{k}}^3(\mathbf{r}).$$

$${\epsilon }_{\text{k}}^1(\mathbf{r})=-\frac{1}{4}{\nabla }^{2}n(\mathbf{r})+\frac{1}{8}\frac{{\left(\nabla n(\mathbf{r})\right)}^2}{n(\mathbf{r})}.$$

$${\epsilon }_{\text{k}}^2(\mathbf{r})=\frac{1}{8}\left(\frac{{\left[\nabla (\rho (\mathbf{r})/n(\mathbf{r}))\right]}^2}{\rho (\mathbf{r})/n(\mathbf{r})}+\frac{{\left[\nabla (S(\mathbf{r})/n(\mathbf{r}))\right]}^2}{S(\mathbf{r})/n(\mathbf{r})}\right)n(\mathbf{r}).$$

$${\epsilon }_{\text{k}}^3\left(\mathbf{r}\right)$$

$$+\frac{1}{4}\sum _{k=1}^3\nabla (S_k(\mathbf{r})/S(\mathbf{r}))\nabla S(\mathbf{r}).$$

$${\epsilon }_{\text{k}}^{\text{LPS}}(\mathbf{r})$$

$${\epsilon }_{\text{k}}^{\text{2D}}(\mathbf{r})$$

$$E_{\text{e-n}}$$

$$E_{\text{e-e}}$$

$$E_{\text{kin}}$$

$$+\frac{1}{8}\int \left(\frac{{\left[\nabla \rho (\mathbf{r})\right]}^2}{\rho (\mathbf{r})}+\sum _{k=1}^3\frac{{\left[\nabla S_k(\mathbf{r})\right]}^2}{S_k(\mathbf{r})}\right)\text{d}\mathbf{r}.$$

$$T_{vw}=\frac{1}{8}\int \frac{(\nabla n{)}^2}{n},$$

$${\rho }^{\text{nicefrac}12}(x)=A\sin kx,S_1^{\text{nicefrac}12}(x)=S^{\text{nicefrac}12}(x)=A\cos kx.$$

$${\epsilon }_{\text{k}}^1(\mathbf{r})$$

$$H\Psi$$

$$=\underset{s_0}{\underbrace{\left(-\frac{1}{2}{\nabla }^2\right){\rho }^{\text{nicefrac}12}+v_0{\rho }^{\text{nicefrac}12}-\sum _l{v}_l{S}_l^{\text{nicefrac}12}}}+\sum _k{\mathbf{ie}}_k\underset{b_k}{\underbrace{\left[\left(-\frac{1}{2}{\nabla }^2\right)S_k^{\text{nicefrac}12}+v_0{S}_k^{\text{nicefrac}12}+v_k{\rho }^{\text{nicefrac}12}-\left(v_i{S}_j^{\text{nicefrac}12}-v_j{S}_i^{\text{nicefrac}12}\right)\right]}}.$$

$$\left(\begin{array}{cccc}{\hat{h}}_0& -v_1& -v_2& -v_3\\ v_1& {\hat{h}}_0& v_3& -v_2\\ v_2& -v_3& {\hat{h}}_0& v_1\\ v_3& v_2& -v_1& {\hat{h}}_0\end{array}\right)\left(\begin{array}{c}{\rho }^{\text{nicefrac}12}\\ S_1^{\text{nicefrac}12}\\ S_2^{\text{nicefrac}12}\\ S_3^{\text{nicefrac}12}\end{array}\right)=\mu \left(\begin{array}{c}{\rho }^{\text{nicefrac}12}\\ S_1^{\text{nicefrac}12}\\ S_2^{\text{nicefrac}12}\\ S_3^{\text{nicefrac}12}\end{array}\right).$$

Full text

Orbital-free density functional theory has the potential to vastly increase the computational efficiency of electronic structure calculations Levy (1979); Levy, Perdew, and Sahni (1984); Pearson, Smargiassi, and Madden (1993); Wesolowski and Wang (2013); Lehtomäki et al. (2014); Karasiev and Trickey (2015). However, a general kinetic energy functional has remained elusive Garcia-Aldea and Alvarellos (2008); Huang and Carter (2010); Shin and Carter (2014); Mi, Genova, and Pavanello (2018); Constantin, Fabiano, and Della Sala (2018a, b). It has even been shown that a perturbative approach to the construction of such a functional would diverge Seidl, Perdew, and Kurth (2000). This has been seen, by more traditional DFT theorists, as an indication that the whole endeavour is ultimately doomed.

While approximations to the kinetic energy functionals exist for specific systems, a universal functional does so far not exist. However, given the large amount of work and ingenuity spent on this research over the last thirty years, it could be that the lack of decisive progress is due not to the intractability of nature at this scale, but due to restrictions in the very foundations on which the theoretical models have been constructed.

We have, in a number of publications (see Pope and Hofer (2017, 2019) and references therein), advocated a two-density model, where the general properties of a many-electron system are represented by mass density $\rho$ and spin density $S$ components. The spin density components in this case are associated with a chiral - and complex - vector field. These two components combined yield the conventional DFT density $n({\bf r})$:

[TABLE]

It has been shown that such a two-density model also complies with theorems equal to the Hohenberg-Kohn theorems for a single density Hohenberg and Kohn (1964); Pope and Hofer (2019). Moreover, it allowed us to define an effective wavefunction combining the roots of these two densities in the following way:

[TABLE]

Here, ${\bf e}_{S}$ a unit vector in the direction of the spin vector field and ${\bf i}{\bf e}_{S}$ is a unit bivector in the plane perpendicular to the spin vector. Such an effective wavefunction will be subject to an altered Hamiltonian, essentially due to the multivector properties of the effective wavefunction Doran and Lasenby (2003). The modified Hamiltonian is then given by (in atomic units):

[TABLE]

Here, ${\bf e}_{b}({\bf r})$ is a unit vector in the direction associated with the potential $v_{b}({\bf r})$, which we termed the bivector potential due to its geometric properties. The approach yields a very compact formulation for the many electron problem under the assumption that ${\bf e}_{S}({\bf r})\approx{\bf e}_{S}$, in other words a constant vector. Then the vector ${\bf e}_{b}({\bf r})={\bf e}_{b}={\bf e}_{S}$ and the general problem is determined by the matrix equation (we omit the explicit dependency on the position in the following):

[TABLE]

where $\mu$ is the chemical potential as in the standard formulation. It can be seen that the bivector potential governs the interaction between the two density components. In case it vanishes the formulation simplifies to the single electron case and a conventional Schrödinger equation for single electrons. The bivector potential can be explicitly derived from the gradients of the two density components. It is Pope and Hofer (2019):

[TABLE]

However, it should be noted that the theoretical framework described by the last two equations is a simplification, as it assumes the unit vector of the spin density fields to be constant. In the following, we shall describe model simulations with this simplified model. The rest of the paper is then devoted to removing the simplification. The exact formulation of the problem, as will be shown, also leads to an exact kinetic energy functional.

We have implemented the simplified model in the self-consistency cycle of the DFT package CASTEP Clark et al. (2005). Similar to the auxiliary functional approach Hasnip and Probert (2015), the model takes the Kohn-Sham Kohn and Sham (1965) density at each stage of the self-consistency cycle, independently minimises it according to the two-density many-body equations and returns a density to the parent Kohn-Sham cycle. Convergence occurs when consecutive densities returned by the two-density model are equivalent and thus the results of the two models are consistent. Since correlation is not yet present in our model, we employ a Local-Density-Approximation (LDA) correlation term Perdew and Zunger (1981).

Figure 1 shows the result of model simulations for eighteen atomic systems. The figure shows the absolute values of the energy differences obtained with standard Kohn-Sham DFT methods and with the two-density model. The differences are in all cases very small. We can thus say that the simplified model already agrees with the total energy values obtained with standard DFT methods. The agreement is very good not only for the atomic systems shown here, but also for the interatomic distances in dimers, and for the bandstructure of semiconductors which is very sensitive to variations of the electron density. These simulations are shown elsewhere Pope and Hofer (2019).

It should be noted that the two-density approach no longer contains exchange functionals as in standard DFT. Within our approach, exchange energy is an emergent property, which is due to the interaction between the two densities. So it seems that two very different approaches to a general many-electron system lead to the same physical properties. Moreover, the agreement, in case it persists for more complicated systems, raises questions about the level of approximation present in standard DFT methods. Because within the two-density model the results were derived under the assumption that the vector fields for both, spin-densities and bivector potentials, are independent of location and of constant direction. In general, this will not be the case.

It is difficult, at this point in the development of the two-density method, to gauge the importance of this finding. Generally, it had been assumed within the DFT community, that only the exchange-correlation functionals are subject to approximations. If we knew their exact form, then we would have an in principle exact method of simulating many-electron systems, was the general opinion so far. This result, however, suggests that there is another level of approximation involved, which does not depend on the exchange-correlation functionals at all. And this level concerns the change of local physical properties. If the direction of the vector associated with the spin density, which is one of the two density components, changes faster than a certain threshold, then standard DFT methods will no longer be sufficient to describe the dynamics of the system. In the following we shall derive the exact functional and show the approximations involved in the two-density scalar kinetic energy, which is equivalent to the standard DFT formalism.

In general both vectors, ${\bf e}_{S}$ and ${\bf e}_{b}$ will depend on the position so that ${\bf e}_{S}={\bf e}_{S}({\bf r})$ and ${\bf e}_{b}={\bf e}_{b}({\bf r})$. This means that both, the roots of the spin density and the bivector potential will consist of three individual components. In a cartesian basis we may write the term associated with the spin density as:

[TABLE]

Correspondingly, the term with the bivector potential will, in the general case, be:

[TABLE]

The general formulation is quite complicated, and the detailed analytical calculation is omitted for brevity. It can be found in the supplemental material at § A.1. The ansatz leads to a set of four differential equations. The equations in matrix notation are given by:

[TABLE]

The components $v_{k}$ of the bivector potential are described by the following expression:

[TABLE]

where the differential operator $\hat{\bf J}_{k}$ is given by the expression:

[TABLE]

and the mixing parameters are given by,

[TABLE]

Here, the indices are defined by the bivector notation, ${\bf i}{\bf e}_{k}={\bf e}_{i}{\bf e}_{j}$.

The formulation leads to a very straightforward evaluation of individual energy components of the system. The total energy density is given by,

[TABLE]

where the same Hamiltonian operator is used on each of the four density components (details may be found in § A.2). This operator consists of a kinetic energy operator and a potential associated with the electron-electron and electron-nuclear Coulomb interactions. This leads to three constituent energy densities, the electron-nuclear energy density,

[TABLE]

the electron-electron energy density,

[TABLE]

and the kinetic energy density,

[TABLE]

This form of the energy density has some similarities to previous forumlations, which become clearer when the equation is reformulated. We define the kinetic energy density as a sum of three distinct terms,

[TABLE]

Here the first term requires only the total charge density, $n$, and no further information,

[TABLE]

This is equivalent to the single density formulation presented by Levi et al (LPS) Levy, Perdew, and Sahni (1984).

The second term contains the mass and spin density scalar quantities, $\rho$ and $S$,

[TABLE]

This is equivalent to the kinetic energy density correction presented in the scalar two-density model (2D) Pope and Hofer (2019).

The final term includes the vector information for the spin density,

[TABLE]

This final term completes the kinetic energy density and reveals the approximations that have been made by the previous approaches. Namely, the LPS and 2D approaches can be written respectively,

[TABLE]

A functional for each energy component is given by integrating the energy densities,

[TABLE]

It is noteworthy that the only energy term, which contains the individual components of the total density, is the kinetic energy. Moreover, it contains these components in a form which is very familiar in orbital-free DFT. The original von Weizsäcker correction to the kinetic energy functional has the form Von Weizsacker (1935):

[TABLE]

where $n$ is the electron (mass or charge) density. It can be seen that the two-density model leads to similar terms. This suggests that the correction has a more general significance than previously assumed. Here, it is not just a term introduced from a single-electron case, that is the solution of the hydrogen problem, but a term which is always present in a many-electron system. This derivation also indicates that it would be largely impossible to capture this behaviour of a many-electron system with a single density only, as has been attempted in the past. An estimate of the exact functional also shows that its kernel is at most linear with the density. The functional for free electrons, by contrast, would involve the DFT density to the power of 5/3. From this perspective it is quite clear that a model starting from free electrons will lead to instabilities as soon as the density increases beyond a certain level.

We note that this kinetic energy functional is the kinetic energy functional of interacting electrons, not non-interacting electrons as in the Kohn-Sham formalism. Furthermore, it is easy to establish that the two-density framework extends into the regime of single free electrons, in contrast to single-density DFT, where the kinetic energy of free electrons cannot be derived from the density. Here, the density roots of free electrons are (we set the spin vector equal to the $x$-direction):

[TABLE]

Then the components of the kinetic energy density read:

[TABLE]

and the kinetic energy functional yields $E_{\text{kin}}=k^{2}/2$, because $A^{2}=1/V$. Notable, the energy derived from a single-density perspective, $\varepsilon_{\text{k}}^{1}$, is zero. From a physical perspective this means that there is continuity in the new framework going from single non-interacting electrons to a more complicated system where electrons interact. The change from the former to the latter can then be as incremental as one would wish. This should make the analysis of events at the local level much simpler in the future.

Acknowledgements

The authors acknowledge EPSRC funding for the UKCP consortium (Grant No. EP/K013610/1). This work was also supported by the North East Centre for Energy Materials (NECEM). WH acknowledges support from the University of Chinese Academy of Sciences.

Appendix A Supporting Material

A.1 Bivector Potential

The Schrödinger equation in the extended electron model is given,

[TABLE]

Here, the indices $i$ and $j$ are defined by bivector notation, such that ${\bf ie}_{k}={\bf e}_{i}{\bf e}_{j}$. This yields four equations, the scalar and the three bivector parts, which in matrix form are given,

[TABLE]

Multiplying by the dual, we find,

[TABLE]

The scaler part may be simplified,

[TABLE]

We note that one can show,

[TABLE]

The bivector part may also be simplified,

[TABLE]

Here we introduce the notation,

[TABLE]

So the bivector term can be written,

[TABLE]

Since we know the total energy must be a scalar object, the bivector terms must be zero,

[TABLE]

where, we define,

[TABLE]

so that,

[TABLE]

substituting the expresion for $v_{j}(v_{i})$ into the equation for $v_{i}(v_{j})$, we find

[TABLE]

and subtituting these expressions into the equation for $v_{k}$ yeilds,

[TABLE]

where

[TABLE]

A.2 Kinetic Energy Density

The total energy density is given by Eq. 12,

[TABLE]

where $v_{0}$ contains the electron-electron and electron-nuclear potentials. The kinetic energy density is thus given,

[TABLE]

Firstly, let $S_{k}=\varphi_{k}S$,

[TABLE]

Next, let $\rho=\chi n$ and $S=\left(1-\chi\right)n$,

[TABLE]

where

[TABLE]