Exact orbital-free kinetic energy functional for general many-electron systems

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.

Graphical abstract

Keywords

condensed matter / density functional theory (DFT) / extended electrons

Orbital-free density functional theory has the potential to vastly increase the computational efficiency of electronic structure calculations [1-6]. However, a general kinetic energy functional has remained elusive [7-12]. It has even been shown that a perturbative approach to the construction of such a functional would diverge [13]. 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 [14, 15] and references therein), advocated a two-density model, where the general properties of a many-electron system are represented by mass density ρ 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 (r):

\begin{array}{}\text{(1)}& n(r)=\rho (r)+S(r)\mathrm{.}\end{array}

It has been shown that such a two-density model also complies with theorems equal to the Hohenberg–Kohn theorems for a single density [15, 16]. Moreover, it allowed us to define an effective wavefunction combining the roots of these two densities in the following way:

\begin{array}{}\text{(2)}& \Psi (r)={\rho }^{\frac{1}{2}}(r)+\text{i}{e}_S(r)S^{\frac{1}{2}}(r)\mathrm{.}\end{array}

Here, e_S a unit vector in the direction of the spin vector field and ie_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 [17]. The modified Hamiltonian is then given by (in atomic units):

\begin{array}{}\text{(3)}& \begin{array}{ll}H(r)\hfill & =-\frac{1}{2}{\nabla }^2+v_0(r)+\text{i}{e}_b(r)v_b(r)\hfill \\ \hfill & =:{\widehat{h}}_0(r)+\text{i}{e}_b(r)v_b(r)\mathrm{.}\hfill \end{array}\end{array}

Here, e_b(r) is a unit vector in the direction associated with the potential ve_b(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 e_S(r)≈e_S, in other words a constant vector. Then the vector e_b(r)=e_b=e_S and the general problem is determined by the matrix equation (we omit the explicit dependency on the position in the following):

\begin{array}{}\text{(4)}& \left(\begin{array}{cc}{\widehat{h}}_0& -v_b\\ v_b& {\widehat{h}}_0\end{array}\right)\left(\begin{array}{c}{\rho }^{\frac{1}{2}}\\ S^{\frac{1}{2}}\end{array}\right)=\mu \left(\begin{array}{c}{\rho }^{\frac{1}{2}}\\ S^{\frac{1}{2}}\end{array}\right)\mathrm{,}\end{array}

where μ 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 [15]

\begin{array}{}\text{(5)}& v_b=\frac{1}{\rho +S}\left[S^{\frac{1}{2}}\left(-\frac{{\nabla }^2}{2}\right){\rho }^{\frac{1}{2}}-{\rho }^{\frac{1}{2}}\left(-\frac{{\nabla }^2}{2}\right)S^{\frac{1}{2}}\right]\mathrm{.}\end{array}

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 [18]. Similar to the auxiliary functional approach [19], the model takes the Kohn-Sham [20] 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 [21].

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 [15].

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, e_S and e_b will depend on the position so that e_S=e_S(r) and e_b=e_b(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

\begin{array}{}\text{(6)}& \text{i}{e}_S(r)S^{\frac{1}{2}}(r)={\displaystyle \sum _{j=1}^3}\text{i}{e}_j{S}_j^{\frac{1}{2}}(r)\mathrm{.}\end{array}

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

\begin{array}{}\text{(7)}& \text{i}{e}_b(r)v_b(r)={\displaystyle \sum _{j=1}^3}\text{i}{e}_j{v}_j(r)\mathrm{.}\end{array}

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

\begin{array}{}\text{(8)}& \left(\begin{array}{cccc}{\widehat{h}}_0& -v_1& -v_2& -v_3\\ v_1& {\widehat{h}}_0& v_3& -v_2\\ v_2& -v_3& {\widehat{h}}_0& v_1\\ v_3& v_2& -v_1& {\widehat{h}}_0\end{array}\right)\left(\begin{array}{c}{\rho }^{\frac{1}{2}}\\ S_1^{\frac{1}{2}}\\ S_2^{\frac{1}{2}}\\ S_3^{\frac{1}{2}}\end{array}\right)=\mu \left(\begin{array}{c}{\rho }^{\frac{1}{2}}\\ S_1^{\frac{1}{2}}\\ S_2^{\frac{1}{2}}\\ S_3^{\frac{1}{2}}\end{array}\right)\mathrm{.}\end{array}

The components v_k of the bivector potential are described by the following expression:

\begin{array}{}\text{(9)}& v_k=\frac{1}{4}\frac{{\widehat{J}}_k-B_k^j{\widehat{J}}_i-B_k^i{\widehat{J}}_j}{1-B_k^j{\chi }_i^{+}-B_k^i{\chi }_j^{-}}\mathrm{,}\end{array}

where the differential operator {\widehat{J}}_k is given by the expression:

\begin{array}{}\text{(10)}& \begin{array}{ll}{\widehat{J}}_k=\hfill & {\left(\frac{\rho -S}{2}+S_k\right)}^{-1}\hfill \\ \hfill & \times \nabla \left({\rho }^{\frac{1}{2}}\nabla S_k^{\frac{1}{2}}-S_k^{\frac{1}{2}}\nabla {\rho }^{\frac{1}{2}}+S_i^{\frac{1}{2}}\nabla S_j^{\frac{1}{2}}-S_j^{\frac{1}{2}}\nabla S_j^{\frac{1}{2}}\right)\mathrm{,}\hfill \end{array}\end{array}

and the mixing parameters are given by

\begin{array}{}\text{(11)}& \begin{array}{l}{B}_k^j=\frac{{\chi }_j^{-}-{\chi }_i^{+}{\chi }_k^{+}}{1-{\chi }_k^{-}{\chi }_k^{+}}\mathrm{,}{B}_k^i=\frac{{\chi }_i^{+}-{\chi }_j^{-}{\chi }_k^{-}}{1-{\chi }_k^{-}{\chi }_k^{+}}\mathrm{,}\\ {\chi }_k^{\pm }=\left(S_i^{\frac{1}{2}}{S}_j^{\frac{1}{2}}\pm {\rho }^{\frac{1}{2}}{S}_k^{\frac{1}{2}}\right){\left(\frac{\rho -S}{2}+S_k\right)}^{-1}\mathrm{.}\end{array}\end{array}

Here, the indices are defined by the bivector notation, ie_k=e_ie_j.

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

\begin{array}{}\text{(12)}& {ϵ}_{\text{tot}}={\rho }^{\frac{1}{2}}{\widehat{h}}_0{\rho }^{\frac{1}{2}}+{\displaystyle \sum _k}{S}_k^{\frac{1}{2}}{\widehat{h}}_0{S}_k^{\frac{1}{2}}\mathrm{,}\end{array}

where the same Hamiltonian operator is used on each of the four density components (details may be found in Appendix 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,

\begin{array}{}\text{(13)}& {\epsilon }_{\text{e}-\text{n}}(r)=v_{\text{e}-\text{n}}(r)n(r)\mathrm{,}\end{array}

the electron–electron energy density,

\begin{array}{}\text{(14)}& {\epsilon }_{\text{e}-\text{e}}(r)=v_{\text{e}-\text{e}}(r)n(r)\mathrm{,}\end{array}

and the kinetic energy density,

\begin{array}{}\text{(15)}& {\epsilon }_k(r)=-\frac{1}{4}{\nabla }^{2}n(r)+\frac{1}{8}\left(\frac{{\left[\nabla \rho (r)\right]}^2}{\rho (r)}+{\displaystyle \sum _{k=1}^3}\frac{{\left[\nabla S_k(r)\right]}^2}{S_k(r)}\right)\mathrm{.}\end{array}

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,

\begin{array}{}\text{(16)}& {\epsilon }_k(r)={\epsilon }_k^1(r)+{\epsilon }_k^2(r)+{\epsilon }_k^3(r)\mathrm{.}\end{array}

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

\begin{array}{}\text{(17)}& {\epsilon }_k^1(r)=-\frac{1}{4}{\nabla }^{2}n(r)+\frac{1}{8}\frac{{[\nabla n(r)]}^2}{n(r)}\mathrm{.}\end{array}

This is equivalent to the single density formulation presented by Levi et al. (LPS) [2].

The second term contains the mass and spin density scalar quantities, ρ and S,

\begin{array}{}\text{(18)}& {\epsilon }_k^2(r)=\frac{1}{8}\left(\frac{{\left[\nabla (\rho (r)/n(r))\right]}^2}{\rho (r)/n(r)}+\frac{{\left[\nabla (S(r)/n(r))\right]}^2}{S(r)/n(r)}\right)n(r)\mathrm{.}\end{array}

This is equivalent to the kinetic energy density correction presented in the scalar two-density (2D) model [15].

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

\begin{array}{}\text{(19)}& \begin{array}{ll}{\epsilon }_k^3\left(r\right)=\hfill & \frac{1}{8}{\displaystyle \sum _{k=1}^3}\frac{{\left[\nabla (S_k(r)/S(r))\right]}^2}{S_k(r)/S(r)}\hfill \\ \hfill & +\frac{1}{4}{\displaystyle \sum _{k=1}^3}\nabla (S_k(r)/S(r))\nabla S(r)\mathrm{.}\hfill \end{array}\end{array}

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,

\begin{array}{}\text{(20)}& \begin{array}{ll}{\epsilon }_k^{\text{LPS}}(r)\hfill & ={\epsilon }_k^1(r)\mathrm{,}\hfill \\ {\epsilon }_k^{\text{2D}}(r)\hfill & ={\epsilon }_k^1(r)+{\epsilon }_k^2(r)\mathrm{.}\hfill \end{array}\end{array}

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

\begin{array}{}\text{(21)}& \begin{array}{lll}{E}_{\text{e}-\text{n}}\hfill & =\hfill & {\displaystyle \int }{v}_{\text{e}-\text{n}}(r)n(r)\text{d}r\mathrm{,}\hfill \\ E_{\text{e}-\text{e}}\hfill & =\hfill & {\displaystyle \int }{v}_{\text{e}-\text{e}}(r)n(r)\text{d}r\mathrm{,}\hfill \\ E_{\text{kin}}\hfill & \text{=}\hfill & -\frac{1}{4}{\displaystyle \int }{\nabla }^{2}n(r)\text{d}r\hfill \\ \hfill & \hfill & +\frac{1}{8}{\displaystyle \int }\left(\frac{{\left[\nabla \rho (r)\right]}^2}{\rho (r)}+{\displaystyle \sum _{k=1}^3}\frac{{\left[\nabla S_k(r)\right]}^2}{S_k(r)}\right)\text{d}r\mathrm{.}\hfill \end{array}\end{array}

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 [22]:

\begin{array}{}\text{(22)}& T_{vw}=\frac{1}{8}{\displaystyle \int }\frac{{(\nabla n)}^2}{n}\mathrm{,}\end{array}

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):

\begin{array}{}\text{(23)}& {\rho }^{\frac{1}{2}}(x)=A\sin kx\mathrm{,}{S}_1^{\frac{1}{2}}(x)=S^{\frac{1}{2}}(x)=A\cos kx\mathrm{.}\end{array}

Then the components of the kinetic energy density read

\begin{array}{}\text{(24)}& {\epsilon }_k^1(r)=\mathrm{0,}{\epsilon }_k^2(r)=\frac{1}{2}{A}^2{k}^2\mathrm{,}{\epsilon }_k^3(r)=0.\end{array}

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, {\epsilon }_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.

References

[1] M. Levy, Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem, Proc. Natl. Acad. Sci. USA 76(12), 6062 (1979)

[2] M. Levy, J. P. Perdew, and V. Sahni, Exact differential equation for the density and ionization energy of a many-particle system, Phys. Rev. A 30(5), 2745 (1984)

[3] M. Pearson, E. Smargiassi, and P. Madden, Ab initio molecular dynamics with an orbital-free density functional, J. Phys.: Condens. Matter 5(19), 3221 (1993)

[4] T. A. Wesolowski and Y. A. Wang, Recent Progress in Orbital Free Density Functional Theory, Vol. 6, World Scientific, 2013

[5] J. Lehtomäki, I. Makkonen, M. A. Caro, A. Harju, and O. Lopez Acevedo, Orbital-free density functional theory implementation with the projector augmented-wave method, J. Chem. Phys. 141(23), 234102 (2014)

[6] V. V. Karasiev and S. B. Trickey, Frank discussion of the status of ground-state orbital-free DFT, in: Advances in Quantum Chemistry, Vol. 71, Elsevier, 2015, pp 221–245

[7] D. García-Aldea and J. Alvarellos, Approach to kinetic energy density functionals: Nonlocal terms with the structure of the von Weizsäcker functional, Phys. Rev. A 77(2), 022502 (2008)

[8] C. Huang and E. A. Carter, Nonlocal orbital-free kinetic energy density functional for semiconductors, Phys. Rev. B 81(4), 045206 (2010)

[9] I. Shin and E. A. Carter, Enhanced von Weizsäcker Wang–Govind–Carter kinetic energy density functional for semiconductors, J. Chem. Phys. 140, 18A531 (2014)

[10] W. Mi, A. Genova, and M. Pavanello, Nonlocal kinetic energy functionals by functional integration, J. Chem. Phys. 148(18), 184107 (2018)

[11] L. A. Constantin, E. Fabiano, and F. Della Sala, Semilocal Pauli–Gaussian kinetic functionals for orbital-free density functional theory calculations of solids, J. Phys. Chem. Lett. 9(15), 4385 (2018)

[12] L. A. Constantin, E. Fabiano, and F. Della Sala, Nonlocal kinetic energy functional from the jellium-with-gap model: Applications to orbital-free density functional theory, Phys. Rev. B 97(20), 205137 (2018)

[13] M. Seidl, J. P. Perdew, and S. Kurth, Simulation of all-order density-functional perturbation theory, using the second order and the strong-correlation limit, Phys. Rev. Lett. 84(22), 5070 (2000)

[14] T. Pope and W. Hofer, Spin in the extended electron model, Front. Phys. 12(3), 128503 (2017)

[15] T. Pope and W. Hofer, A two-density approach to the general many-body problem and a proof of principle for small atoms and molecules, Front. Phys. 14(2), 23604 (2019)

[16] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136(3B), B864 (1964)

[17] C. Doran and A. Lasenby, Geometric Algebra for Physicists, Cambridge University Press, 2003

[18] S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. Probert, K. Refson, and M. C. Payne, First principles methods using CASTEP, Z. Kristallogr. Cryst. Mater. 220(5/6), 567 (2005)

[19] P. Hasnip and M. Probert, Auxiliary density functionals: a new class of methods for efficient, stable density functional theory calculations, arXiv: 1503.01420 (2015)

[20] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140(4A), A1133 (1965)

[21] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23(10), 5048 (1981)

[22] C. Von Weizsacker, On the theory of nuclear masses, Z. Phys. 96, 431 (1935)

Rights & permissions

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature 2020