Orbital-Free Density Functional Theory: Kinetic Potentials and Ab Initio Local
Pseudopotentials
Jeng-Da Chai∗1 and John D. Weeks1, 2
Institute for Physical Science and Technology,
and Department of Chemistry and Biochemistry,
University of Maryland, College Park, MD 20742
(Dated: November 1, 2018)
In the density functional (DF) theory of Kohn and Sham, the kinetic energy of the ground state
of a system of noninteracting electrons in a general external field is calculated using a set of orbitals.
Orbital free methods attempt to calculate this directly from the electron density by approximating
the universal but unknown kinetic energy density functional. However simple local approximations
are inaccurate and it has proved very difficult to devise generally accurate nonlocal approximations.
We focus instead on the kinetic potential, the functional derivative of the kinetic energy DF, which
appears in the Euler equation for the electron density. We argue that the kinetic potential is more
local and more amenable to simple physically motivated approximations in many relevant cases, and
describe two pathways by which the value of the kinetic energy can be efficiently calculated. We
propose two nonlocal orbital free kinetic potentials that reduce to known exact forms for both slowly
varying and rapidly varying perturbations and also reproduce exact results for the linear response of
the density of the homogeneous system to small perturbations. A simple and systematic approach
for generating accurate and weak ab initio local pseudopotentials which produce a smooth slowly
varying valence component of the electron density is proposed for use in orbital free DF calculations
of molecules and solids. The use of these local pseudopotentials further minimizes the possible errors
from the kinetic potentials. Our theory yields results for the total energies and ionization energies
of atoms, and for the shell structure in the atomic radial density profiles that are in very good
agreement with calculations using the full Kohn-Sham theory.
I. INTRODUCTION
Density-functional theory (DFT) has become one of
the most powerful tools for investigating the electronic
structure of large complex systems. In principle, as
shown by Hohenberg and Kohn [1], the exact ground
state energy of a system of N electrons can be formally
written as a functional E[ρ] of only the electron density
ρ(r), a function of three variables, and the external field
Vext(r). Determining the energy and other ground state
properties from such an approach could dramatically
reduce the computational cost for large systems when
compared with traditional quantum chemistry methods,
which deal with wavefunctions involving coordinates of
all N electrons [2, 3].
Kohn and Sham (KS) [4, 5] showed that E[ρ] can be
usefully partitioned into the following set of terms:
E[ρ] = Ts[ρ] + EH [ρ] + Exc[ρ] +
ρ(r)Vext(r)dr. (1)
Here Ts[ρ] is the noninteracting kinetic energy density
functional (KEDF), which gives the kinetic energy of
a model system of N noninteracting electrons in a self-
consistent field chosen so that the ground state density
∗E-mail: jdchai@berkeley.edu. Present address: Molecular
Foundry, Materials Sciences Division, Lawrence Berkeley National
Laboratory, and Department of Chemistry, University of Califor-
nia, Berkeley, California 94720
equals ρ(r),
EH [ρ] ≡
ρ(r)ρ(r′)
|r− r′|
drdr′ (2)
is the classical electron-electron potential energy (Hartree
energy) and Exc[ρ] is the exchange-correlation energy (in-
cluding the difference between the interacting and non-
interacting kinetic energy and the difference between
the quantum and classical electron-electron potential en-
ergy). The last term on the right of Eq. (1) is the only
term that depends explicitly on the external potential
Vext(r). Atomic units are used throughout the paper.
If all these functionals were known, then the density
ρ(r) could be obtained from the variational principle (Eu-
ler equation) associated with minimizing Eq. (1):
µ = VTs(r; [ρ]) + Veff (r; [ρ]), (3)
and the total energy of the inhomogeneous system could
then be determined from the energy functional E[ρ]. All
other physical quantities related to the ground-state den-
sity could also be computed. Here µ is the chemical po-
tential (the Lagrange multiplier associated with the nor-
malization condition
ρ(r)dr = N), and Veff (r; [ρ]) is
an effective one-body potential defined by
Veff (r; [ρ]) ≡
δρ(r)
EH [ρ] + Exc[ρ] +
ρ(r)Vext(r)dr
= VH(r; [ρ]) + Vxc(r; [ρ]) + Vext(r), (4)
where
VH(r; [ρ]) ≡ δEH [ρ]/δρ(r) =
ρ(r′)
|r− r′|
dr′ (5)
http://arxiv.org/abs/0704.1878v1
is the Hartree potential, and Vxc(r; [ρ]) ≡ δExc[ρ]/δρ(r) is
the exchange-correlation potential. Similarly we interpret
VTs(r; [ρ]) ≡ δTs[ρ]/δρ(r) (6)
as the kinetic potential (KP) arising from the KEDF [6].
Further progress requires an accurate determination
of the noninteracting kinetic energy, whose magnitude is
much larger than the exchange-correlation energy. The
initial development of DFT as a practical computational
method was made possible by KS’s realization that the
numerical value of the noninteracting kinetic energy can
be exactly calculated, not directly from the density itself
using Ts[ρ], but by introducing a set of N one-electron
wave functions (orbitals) satisfying the N coupled KS
equations that describe the model system [4, 5].
Research could then focus on determining the remain-
ing small term Exc[ρ]. Here even local density approx-
imations have often proved useful. Through the efforts
of many workers we now have generally accurate expres-
sions for Exc[ρ]. Using these along with the KS orbitals
to calculate the kinetic energy, one can accurately cal-
culate both the total energy E[ρ] and the ground-state
density ρ(r) for a wide variety of systems.
However, the use of the KS orbitals usually generates
a relatively expensive O(N3) scaling of computational
cost with the number of electrons. While this scaling is
much better than that of most standard methods that
include correlation energy, calculations for large systems
remain problematic. This remaining bottleneck could be
removed if there were an accurate treatment of the kinetic
energy in terms of the electron density only [2, 3, 7].
To that end there has been considerable effort in-
vested in developing “orbital-free”density functional the-
ory (OF-DFT) by making direct approximations for Ts[ρ]
[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. While
earlier simple local density approximations for Ts[ρ] like
those used in the Thomas-Fermi (TF) model [21] are very
inaccurate, there have been two main advances in recent
work that offer prospects for significant improvements.
The first is the introduction of nonlocal KEDFs that
reproduce known exact results for very slowly varying or
very rapidly varying fields and give the exact linear re-
sponse (LR) of the density of the uniform model system
to small perturbations. Similar ideas have been success-
fully applied to classical nonuniform fluids [22]. The sec-
ond advance is to focus not on the total density but on
the smaller and more slowly varying valence electron den-
sity as described by a weak pseudopotential acting only
on the valence electrons. While conventional pseudopo-
tential methods use orbitals, recently developed ab ini-
tio local pseudopotential (AILPS) methods determine the
unique local one-body potential producing a given target
valence density by solving the KS equations inversely, us-
ing the one-to-one mapping between density and poten-
tials in DFT [23]. For OF-DFT with LR-based KEDF’s,
the use of pseudopotentials not only can reduce the com-
putational cost, but also can improve its accuracy, since
the system will be closer to the LR regime where Ts[ρ] is
designed to be accurate. [24, 25]. Indeed, very promising
results using such OF-DFT methods have been obtained
for a variety of nearly free-electron-like metals.
However, existing KEDF’s have not yet achieved chem-
ical accuracy for systems with localized and more rapidly
varying electron densities like molecules or for covalent
or ionic solids. The main problem is that the exact Ts[ρ]
is highly nonlocal, and we have little idea of the func-
tional form of the nonlocality for densities far from the
LR regime. It has proved very difficult to understand
what errors an approximate nonlocal Ts[ρ] will produce
in the density as determined by the Euler equation with
a general Vext(r).
We explore here a different way to attack this basic
problem. The exact Ts[ρ] can be formally obtained from
VTs(r; [ρ]) by functional integration over density changes
in all regions of space [6, 18]. Because of this integration
Ts[ρ] is a more nonlocal functional of the density than
is VTs(r; [ρ]). More detailed arguments arriving at this
same conclusion have been recently presented [26]. Since
most problems in devising accurate approximations for
Ts[ρ] have arisen from the nonlocality, this suggests it
could be worthwhile to try to develop approximations
for the KP VTs(r; [ρ]) itself.
To illustrate this point, Chai and Weeks [7] added a
simple gradient correction to the original local TF KP
for atoms [21], with a coefficient chosen to reproduce
the exact boundary condition of exponential decay of the
electron density far from the nucleus. Though quanti-
tative results were not obtained, the resulting modified
Thomas-Fermi (MTF) model gave energies for atoms and
for closed-shell diatomic molecules that showed notable
improvements when compared to the original TF and re-
lated gradient corrected KEDF models. However, the
local gradient correction used in the MTF KP cannot re-
produce the oscillatory atomic shell structure and it does
not satisfy the exact LR behavior in the homogeneous
limit. It is clear that nonlocality even in the KP must be
taken into account to achieve more accurate results.
We propose here new nonlocal approximations for the
KP using ideas similar to those employed for the nonlocal
KEDFs. These new KPs satisfy the exact LR condition
in the uniform limit, and reproduce known exact limiting
forms of VTs(r; [ρ]) both for very slowly varying and very
rapidly varying perturbations. As will become clear, the
nonlocality in our KP is determined by the requirement
that LR is exactly satisfied, and it is much easier to en-
sure that LR holds for the KP than it is for analogous
KEDF models. We believe this level of nonlocality in the
KP may suffice in many cases when used in conjunction
with AILPS methods to describe slowly varying valence
density components closer to the LR regime.
The remainder of this paper is organized as follow.
Section II will discuss some general pathways connect-
ing Ts[ρ] and VTs(r; [ρ]). Section III will describe limit-
ing forms of the KEDF and KP for slowly varying and
rapidly varying perturbations, and discuss LR theory, an
exact theory for the response of density of the uniform
electron gas to small perturbations. Section IV will de-
velop two nonlocal KPs incorporating both the correct
limiting forms of the exact KP and the exact LR of the
free-electron gas. Section V will compare the numerical
results of the present method for atoms with the KS-DFT
and other KEDFs, both for all-electron calculations and
for valence electrons using the AILPS. We find that the
use of AILPS indeed reduces errors arising from nonlocal-
ity in these approximate KPs or KEDFs, which give very
accurate results for the relatively slowly varying valence
densities. Our conclusions are given in Section VI.
II. PATHWAYS FROM VTs(r; [ρ]) TO Ts[ρ]
If Ts[ρ] is known, VTs(r; [ρ]) can be simply computed
by functional differentiation. However, there is no unique
way of determining Ts[ρ] from a given VTs(r; [ρ]). Many
possible pathways can be used to construct Ts[ρ] by func-
tional integration of VTs(r; [ρ]) [27, 28]. If the exact
VTs(r; [ρ]) is used and the integration is carried out ex-
actly, then all pathways would give the same exact result
for Ts[ρ]. However, when an approximate VTs(r; [ρ]) is
used, different pathways will give different results for the
kinetic energy. But this “thermodynamic inconsistency”
is small if reasonably good approximations are used, since
the integration tends to smooth out local errors that may
exist in the density [7, 28].
More problematic is the fact that most pathways re-
quire additional results for partially coupled systems as
the external field or density perturbation is gradually
turned on, which adds to the computational burden. In
particular, most earlier work has used a “potential en-
ergy pathway”, where the external potential is scaled by
a coupling parameter [7, 27]. The kinetic energy can then
be found by subtracting the potential energy (calculated
from the potential energy density functionals) from the
total energy. However, this pathway is expensive, since
one has to solve the Euler equation (3) for each partially
coupled V λext(r) (with the same µ), to determine the cor-
responding ρλ.
A. Herring’s Pathway
However, Herring showed there is a particular path-
way arising from exact scaling relations between the non-
interacting kinetic energy Ts[ρ] with respect to the coor-
dinate r in ρ(r) where very simple results involving only
the final density can sometimes be found [6, 18]. If the
coordinate r is scaled to αr, the normalized scaled den-
sity is ρα(r) = α
3ρ(αr). It is easy to show that the exact
Ts[ρ] then obeys
α2Ts[ρ] = Ts[ρα]. (7)
For isolated systems, such as atoms and molecules, the
density and its derivatives to all order vanish far from the
nuclei. For such systems, when Eq. (7) is differentiated
with respect to α, and the partial derivative is evaluated
at α = 1, we find the formally exact result
Ts[ρ] =
VTs(r; [ρ])∇ · (rρ(r))dr. (8)
Therefore, once the kinetic potential VTs(r; [ρ]) is
known for some given ρ(r), the numerical value of Ts[ρ]
can then be immediately determined from Eq. (8). Since
there is no need to perform a coupling parameter integra-
tion over the change of density or potential, this scheme
is not only fast, but also numerically reliable. The final
form of Eq. (8) is essentially the viral theorem, and is
directly related to the force on molecules [29, 30, 31, 32].
Note that this simple and exact pathway holds only
for the noninteracting Ts[ρ] [31], which again shows the
virtues of the KS partitioning of the total energy. We will
use Eq. (8) as the basic pathway to determine the numer-
ical value of Ts[ρ] from a given approximate VTs(r; [ρ])
for most calculations in this paper. However, Eq. (8)
does not hold for extended solid-state systems because of
nonvanishing boundary terms, and thus far we have not
found an exact and simple way of including them.
Fortunately, there is another class of computationally
efficient“density pathways”that can be used for extended
systems, as we now show. Density pathways can also
be used for atomic and molecular systems to check the
accuracy of the VTs used, since results using the exact
VTs would be independent of path [27, 28].
B. Density Pathways
The change in the kinetic energy can be formally re-
lated to a coupling parameter integration, where the den-
sity changes from some known value at λ = 0 to the final
density at λ = 1:
Ts[ρ] = Tλ=0 +
drVTs(r; [ρλ])
∂ρλ(r)
In most cases a simple linear density pathway will suffice.
Here the density ρ(r) is linearly scaled by a coupling pa-
rameter λ from some uniform reference density ρ0 natu-
rally chosen to be the uniform electron density N/V in
extended systems:
ρλ(r) = ρ0 + λ[ρ(r) − ρ0]. (10)
Then Eq. (9) becomes
Ts[ρ] = Tλ=0 +
drVTs(r; [ρλ])[ρ(r) − ρ0] (11)
Here Tλ=0 is the kinetic energy of the uniform system,
i.e., the Thomas-Fermi kinetic energy TTF [ρ0]. For ex-
tended systems, where the Herring’s pathway cannot be
used, this density pathway appears to be a good way to
compute T . Other density pathways, like the square-root
pathway introduced by Chen and Weeks [28] to describe
nonuniform hard sphere fluids, can be defined, and have
proved useful in certain applications, but we do not con-
sider them here.
Note from Eq. (10) that ρλ(r) depends only on the fi-
nal density, so evaluation of Eq. (11) is straightforward
and this pathway is computationally efficient. Unlike the
potential energy pathway where the external potential is
scaled, there is no need to solve the Euler equation (3)
for its corresponding external potential V λext(r) at each λ.
However, for isolated systems, where ρ0 = 0, this path-
way is likely to be less accurate than Herring’s pathway,
since it does not automatically satisfy the virial theorem.
III. EXACT LIMITS AND LINEAR RESPONSE
THEORY
Although the exact Ts[ρ] is still unknown, several limit-
ing forms have been discovered for particular density dis-
tributions. These provide important cornerstones that
can be used to construct accurate KEDFs and KPs in
many cases, as will be seen below.
In particular, the Thomas-Fermi (TF) KEDF [21] is
known to be exact for a uniform system:
TTF [ρ] = CF
ρ5/3(r)dr, (12)
where CF =
(3π2)2/3. The TF KEDF TTF [ρ] is de-
rived by local use of uniform free-electron gas model, and
is exact for a system with an infinite number of electrons.
The corresponding expression for the TF KP is
VTF (r; [ρ]) ≡ δTTF [ρ]/δρ(r) =
2/3(r). (13)
This depends only on the local value of ρ2/3(r) and thus
formally is more local than the TF KEDF, whose func-
tional dependence on ρ involves the density at all r. Of
course in this simple case the functional integration of
Eq. (13) can be carried out exactly to yield Eq. (12),
but this cannot be done in general and the nonlocality of
Ts[ρ] has proved problematic.
Results for nonuniform systems are best described in
Fourier space. For a very slowly varying perturbation of
the density, the second-order gradient expansion is exact
[33]. It is easy to see that results correct to second order
at small wavevectors are given by
TTF 1
W [ρ] ≡ TTF [ρ] +
TW [ρ], (14)
where
TW [ρ] ≡
|▽ρ(r)|2
dr (15)
is the von Weizsäcker (W) KEDF [34].
TW [ρ] is exact for a system with one or two electrons,
or where the density can be accurately described by a sin-
gle orbital. Moreover it has been argued [8, 9, 10, 12] that
TW [ρ] gives the correct leading order term for a rapidly
varying perturbation with only high wavevector compo-
nents and that the next order correction is reproduced
TW− 3
TF [ρ] ≡ TW [ρ]−
TTF [ρ]. (16)
The W KP is
VW (r; [ρ]) ≡ δTW [ρ]/δρ(r) =
|▽ρ(r)|2
ρ2(r)
▽2ρ(r)
If we represent the full density by an effective single or-
bital function ψ(r),
ρ(r) = |ψ(r)|2 (18)
then the W KP can be written in a compact form that
will later prove useful:
VW (r; [ρ]) = −
▽2ψ(r)
2ψ(r)
. (19)
Finally, the linear response of the density of a uniform
non-interacting electron gas with density ρ0 to a small
perturbation δV (k) = ǫke
ik·r is exactly known [35],
δρ(k) = χL(q)δV (k). (20)
q ≡ k/2kF (21)
is a dimensionless wavevector, where
kF ≡ (3π
1/3 (22)
is the Fermi wavevector (FWV) and k ≡ |k|. The LR
function χL(q) has the form
χL(q) = −
F−1L (q)
1− q2
1 + q
, (23)
where
FL(q) ≡
1− q2
1 + q
has been called the Lindhard function [10].
It is known that the weak logarithmic singularity at
q = 1 in F−1L (q) is responsible for Friedel oscillations, and
may also be important for the appearance of atomic shell
structure. This singularity further divides the Lindhard
function into two branches in Fourier space: the low-
momentum (q < 1 ) or the low-q (LQ) branch, and the
high-momentum (q > 1) or the high-q (HQ) branch [10].
Figure 1: Linear response functions of a uniform system of
noninteracting Fermions as given by the TF, W, and MTF
(see Ref. [7]) models.
 0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2  2.2  2.4  2.6  2.8  3
exact
The dimensionless response function arising from the
TF KEDF is FTF (q) = 1, and that from the W KEDF is
FW (q) = 3q
2 [36]. Clearly, no linear combination of the
TF and the W KEDFs can reproduce the exact Lindhard
function in Eq. (24). This has the following two limits
[10],
FL(q) =
1 + q2/3 +O(q4) q ≪ 1
3q2 − 3/5 +O(q−2) q ≫ 1
It should be noted that the expansions for both the
low-q and high-q limits are correct to all orders in per-
turbation theory, but valid only in the appropriate limits
in Fourier space. On the other hand, the LR theory is
valid for all wavevectors, but is only accurate for small
perturbations. Therefore, the regime where the response
functions of the two limiting KEDFs deviate from the
exact LR function gives an indication of the range of
wavevectors where the two limiting forms are inaccurate.
As shown in Fig. 1, the response function F−1TF (q) has
no momentum dependence and is only exact at q = 0.
The response function F−1W (q) is exact asymptotically at
high q, and remains fairly accurate for q & 2, but is di-
vergent in the low-q branch, and fails completely for the
nearly uniform electron gas. In contrast, the MTF model
[7] gives a reasonably accurate average description of the
exact response function, especially in the important re-
gion near the singularity at q = 1.
IV. CONSTRUCTION OF NONLOCAL
KINETIC POTENTIALS
A. Kinetic Energy Density Functionals TTFλW [ρ]
and TWλTF [ρ]
Simple linear combinations of the two limiting KEDF’s
in Eqs. (12) and (15), such as the TFλW KEDF [33, 37,
38, 39]
TTFλW [ρ] ≡ TTF [ρ] + λTw[ρ] (26)
and the WλTF KEDF [40, 41, 42, 43]
TWλTF [ρ] ≡ TW [ρ] + λTTF [ρ] (27)
have been widely studied for several decades. The value
of the parameter λ was either determined empirically for
getting good atomic energy or obtained by some semi-
classical arguments.
The advantage of these approaches is the ability to gen-
erate a family of simple KEDF’s easily. It has been shown
empirically that the TF1/5W model can give good val-
ues for atomic energies, but the predicted density profiles
are generally not very accurate, both near and far away
from the nucleus. The TTFλW [ρ] and TWλTF [ρ] function-
als give the correct leading term in the density response
to a slowly-varying perturbation and a rapidly-varying
perturbation respectively, and with particular choices of
λ as in Eqs. (14) and (16) they can reproduce the next
order term. Unfortunately, they then will have an in-
correct leading term in the opposite limit, unless λ = 1.
However it has been shown that TTFW [ρ] with λ = 1
always overestimates the exact Ts[ρ] for various systems
[10]. Finally, none of these functionals can reproduce the
exact response function FL(q) in the homogeneous limit.
Since these models fail to satisfy all the known limits,
and nonlocality in Ts[ρ] is not correctly described, it is
also not surprising that atomic shell structure is missing
in these approaches.
B. Combining TF and W Kinetic Potentials
We argue that it may be more profitable to take advan-
tage of known limiting forms of the KP, rather than the
KEDF, and develop approximations for the more local
VTs(r; [ρ]) directly. Again we can rely on known results
in the linear response regime when the density variations
are not too large.
From Eqs. (14) and (16), the following linear combina-
tions of the TF KP and the W KP in Eqs. (13) and (17)
can reproduce exact results to second order for very small
and very large wavevector perturbations respectively:
VTs(r; [ρ]) ≈
VTF (r; [ρ]) +
VW (r; [ρ]) q ≪ 1
VW (r; [ρ])−
VTF (r; [ρ]) q ≫ 1
Since VTF (r; [ρ]) and VW (r; [ρ]) are the only components
up to second order of the two exact limiting forms of the
KP, we can combine them in analogy to the TFλW and
WλTF models and arrive at generalized KPs.
However, instead of combining them using a fixed pa-
rameter λ, it seems natural to represent them in Fourier
space and allow a wavevector dependence in λ = λ(q) to
connect the limiting forms. The λ(q) can then be chosen
in a very simple way so that the exact LR function is
Figure 2: Weight function f̂(q) for the HQ and LQ KPs.
 0  0.5  1  1.5  2  2.5  3
reproduced for a uniform system with density ρ0. In this
way the LR function bridges the exact limits at large and
small wavevectors, and if the theory is applied to weak
perturbations in the linear response regime for interme-
diate wavevectors we can expect very accurate results.
Here, we derive such generalized KPs based on the KP
for the WλTF model.
C. HQ Kinetic Potential
In analogy to the WλTF model in Eq. (27) we look for
a kinetic potential of the form
V 0HQ(k) = VW (k) + λHQ(q) VTF (k)
= VTF (k) + VW (k) + f̂(q) VTF (k), (29)
where
q = k/2kF (30)
is a dimensionless wavevector normalized by the FWV kF
in Eq. (22) of a uniform reference system with density
ρ0 and f̂(q) = λHQ(q) − 1. The superscript 0 in V
indicates use of a uniform reference system. For a small
perturbation, we can linearize the V 0HQ(k) in Eq. (29).
Requiring that it satisfy LR exactly then determines the
weight function f̂(q) as
f̂(q) = FL(q)− 3q
2 − 1. (31)
See Fig. 2.
We refer to Eq. (29) with Eq. (31) as the HQ KPmodel.
It reproduces the correct high-q limit in Eq. (28) up to the
second order. However, unlike Eq. (16), it also satisfies
the correct low-q limit to leading order and gives exact
results for all q in the linear response regime. Inverse
Fourier transform of Eq. (29) then gives
V 0HQ(r; [ρ], kF ) = VTF (r; [ρ]) + VW (r; [ρ])
f(|r− r′|; kF ) VTF (r
′; [ρ])dr′.
This expression is directly useful for extended systems
where a reasonable ρ0 can be defined. For isolated sys-
tems such as atoms and molecules where the density van-
ishes far from the nuclei, it seems natural to replace kF
in Eq. (32) by the local Fermi wavevector (LFWV)
kF (r) ≡ (3π
2ρ(r))1/3, (33)
though errors may be introduced for rapidly varying den-
sity distributions. Using Eq. (13), this yields the general
form of our proposed HQ kinetic potential:
VHQ(r; [ρ], kF (r)) = VTF (r; [ρ]) + VW (r; [ρ])
f(|r− r′|; kF (r)) ρ
2/3(r′)dr′.
Note that the last term in Eq. (34) is most easily com-
puted in Fourier space as
3(2π)3
f̂(k/2kF (r))ρ
2/3(k)e−ik·rdk, (35)
As can been seen in Eq. (31), the weight function
f̂(k/2kF (r)) is determined analytically. Unlike the LR-
based KEDF approaches, no first-order differential equa-
tion is needed to solve for the weight function in Fourier
space. For atomic systems where the LFWV is used, the
convolution in Eq. (35) must be carried out numerically,
which will lead to a quadratic scaling of the HQ model
(and related LQ model described below) in the numberN
of electrons. For extended systems where one can expand
about the local density ρ0, one can use fast Fourier trans-
forms (FFT’s) for a much more efficient computation of
this integral [44].
D. LQ Kinetic Potential
In analogy to the TFλW model in Eq. (26), we could
similarly generate a KP that is accurate to second order
at low q while still reproducing the leading term at high
q. However this is numerically less useful because the
analogue of Eq. (35) involves the Fourier transform of
VW (r
′; [ρ]), which cannot be simply expressed in terms
of the density. Instead, by empirically taking a properly
chosen component of the density outside the integral we
find that Eq. (34) can be modified to produce a new LQ
KP that is accurate to second order at low q and first
order at high q [45]:
VLQ(r; [ρ], kF (r)) = VTF (r; [ρ]) + VW (r; [ρ])
1/6(r)
f(|r− r′|; kF (r)) ρ
1/2(r′)dr′
Extending these ideas we have constructed a modified
KP that satisfies LR everywhere and is accurate to sec-
ond order at both low and high q [45]. However the
functional form is much more complicated, and little ad-
ditional accuracy is gained from the improved behavior
at very small or very large wavevectors, since all forms
use LR to interpolate for intermediate wavevectors, and
this is where most errors arise in practice. Thus we will
report results here only for the HQ and LQ models.
V. RESULTS FOR ATOMS USING THE HQ
AND LQ MODELS
For completeness and to compare to earlier work, we
first briefly discuss all-electron calculations using the pro-
posed HQ and LQ model KPs and the full atomic poten-
tials. We then describe results using AILPS methods.
These are compared with the KS-DFT, the TFλW mod-
els, and the CAT model introduced by Chacón, Alvarel-
los, and Tarazona [12, 13, 14]. The CAT model is a
LR-based KEDF method, which gives some indication of
shell structure. We employed the latest version, which
uses a nonlocal two-body Fermi wave vector with a pre-
scribed functional form depending on an empirical pa-
rameter β = 1/2 (defined in Eq. (3) of Ref. [14]). This
caused the numerical calculations [46] to be consider-
ably more costly than those of the LQ or HQ models,
which used the local Fermi wave vector as in Eq. (33).
All calculations are spin-restricted and use the local den-
sity approximation (LDA) [47, 48, 49] for the exchange-
correlation functional.
A. All-electron calculations
All-electron calculations consider the density response
to the large and rapidly varying nuclear potential. Since
the system is far from the linear response regime, quanti-
tative results from the HQ and LQ models (or from LR-
based KEDF methods) cannot be expected. However, by
incorporating exact results for very large and very small
wavevectors, these models do correct major deficiencies
of the purely local TF model (which, e.g., predicts an
infinite density at the nucleus!) and even give some qual-
itative indications of atomic shell structure.
The numerical method use the Pauli kinetic potential
V P (r; [ρ]) [50, 51], defined as
V P (r; [ρ]) ≡ VTs(r; [ρ]) − VW (r; [ρ]). (37)
Since VW (r; [ρ]) is the exact KP for a system where the
density can be accurately described by a single orbital,
if V P (r; [ρ]) is omitted, one would essentially obtain the
ground state density of the corresponding Boson system,
where all the electrons are in the same orbital. If we
represent the full density by a single orbital function ψ(r),
so that ρ(r) and VW (r; [ρ]) can be written in the forms of
Eq. (18) and Eq. (19) respectively, we can then combine
V P (r; [ρ]) with the one-body potential Veff (r; [ρ]) in Eq.
(3), and derive a Schrödinger-like equation for the Bose
orbital ψ(r),
▽2 + Veff (r; [ρ]) + V
P (r; [ρ])
ψ(r) = µψ(r). (38)
Table I: Atomic energy E using the KS, LQ, HQ, CAT and
the TFλW models in all-electron calculations. MAPE, the
mean absolute percentage error (relative to the KS method) of
various OF models are given at the bottom of their respective
columns.
KS LQ HQ CAT TF1/5W TFW
He −2.834 −2.565 −2.437 −2.675 −2.911 −1.559
Ne −128.2 −134.3 −126.6 −126.2 −129.5 −86.40
Ar −525.9 −545.9 −512.2 −515.1 −526.2 −375.5
Kr −2750 −2805 −2621 −2712 −2748 −2099
Xe −7229 −7306 −6844 −7141 −7214 −5701
Be −14.45 −14.39 −13.64 −14.11 −14.71 −8.699
Mg −199.1 −207.9 −195.7 −195.2 −200.0 −136.4
C −37.42 −38.97 −36.85 −37.25 −38.41 −24.01
N −54.02 −56.71 −53.59 −53.84 −55.39 −35.33
O −74.47 −78.39 −74.02 −74.08 −76.11 −49.39
Si −288.2 −300.4 −282.5 −282.2 −288.9 −200.5
P −339.9 −354.1 −332.7 −332.8 −340.6 −238.3
S −396.7 −412.8 −387.7 −388.5 −397.3 −279.9
MAPE 4.06% 3.42% 1.82% 1.10% 32.0%
Table II: Electron density at the nucleus ρ(0), using the KS,
LQ, HQ, CAT and the TFλW models in all-electron calcula-
tions. MAPE, the mean absolute percentage error (relative to
the KS method) of various OF models are given at the bottom
of their respective columns.
KS LQ HQ CAT TF1/5W TFW
He 3.525 3.088 2.742 3.600 18.23 0.9515
Ne 614.5 576.6 517.6 613.2 2596 169.6
Ar 3819 3642 3282 3812 1.548 × 104 1093
Be 34.86 30.49 27.17 33.75 158.2 8.952
Mg 1086 1024 920.9 1083 4519 303.0
C 126.0 113.3 101.2 122.8 547.9 33.07
N 203.9 185.6 166.1 200.0 876.7 54.24
O 308.6 284.1 254.6 304.6 1317 83.19
Si 1754 1662 1495 1749 7218 493.9
P 2173 2062 1857 2167 8901 614.5
S 2654 2523 2272 2647 1.083 × 104 753.5
MAPE 7.61% 17.2% 1.14% 331% 72.6%
In other words, Veff (r; [ρ]) + V
P (r; [ρ]) is now the one-
body effective potential for the corresponding Boson sys-
tem with the same electron density. This reduction of
an N -fermion problem to a Boson form is widely imple-
mented in OF-DFT due to its numerical stability and its
easy implementation using existing KS-DFT codes [10].
The associated Pauli potentials for the HQ and LQ
models are immediately obtained by subtraction of the
W KP from Eq. (34) and Eq. (36) respectively. The stan-
dard finite difference method for solving Euler equations
for the TFW models [52] are implemented for the LQ
and HQ models, and the nonlocal terms are evaluated by
Fourier transforms. The choices of radial grids for both of
the real and Fourier space and other detailed numerical
methods are given in Ref. [7]. The kinetic energy for HQ
and LQ models is computed using the Herring pathway
in Eq. (8).
As shown in Table I, the atomic energy calculated by
Figure 3: (Color online) Radial density r2ρ(r) of the Kr atom
using the KS method, the LQ and HQ models, the CAT
model, and the TFλW models (see the inset) with the full
nuclear potential.
 0  0.5  1  1.5  2
 0  0.5  1  1.5  2
TF1/5W
the energy-optimized TF1/5W model is very close to the
KS-DFT, and outperforms all the LR-based models, and
other TFλW models. In Table II, we compare the elec-
tron density at the nucleus ρ(0) for various models. The
TF1/5Wmodel overestimates ρ(0) by about a factor of 4,
while the TFW model underestimates it by about 30%.
The predicted values of ρ(0) for all the LR-based mod-
els are very close to the KS results, and are much better
than the TFλW models.
In Fig. (3), we compare the radial density distribu-
tion r2ρ(r) of the LQ and HQ models to that predicted
by other theories for the Kr atom. Both the TF1/5W
and TFW models predict smooth and structureless ra-
dial density profiles. Using the full Coulombic potential,
all the LR-based models can predict an incipient shell
structure for heavy atoms (Z & 30), and these results
are typical. Since the potential is certainly far beyond
the LR regime, these qualitative results with some sug-
gestion of shell structure are about as good as could be
hoped for. The surprisingly good total energies given in
Table I for the TF1/5W model and the LR-based models
shows that averaged thermodynamic properties are less
sensitive to errors in the KP than is the density profile.
The difference in the results for the LQ, HQ and CAT
models indicates that the LR-based OF theory is being
used outside its range of validity. As shown below, we
gain a significant improvement by using the AILPS to
deal with these difficulties.
B. Ab initio local pseudopotential calculations
As discussed earlier, the use of pseudopotentials in non-
local LR-based OF-DFT can improve the accuracy of the
theory because the weaker pseudopotential is more nearly
in the LR regime, where the theory is designed to be ac-
Figure 4: The smooth target density ρ̃v(r) from Eq. (40), with
parameters given in Table IV for the Si pseudoatom used in
the inverse-KS process, and the valence density ρv(r) pre-
dicted by the LQ and HQ models using the Vps(r) (see Fig.
5) corresponding to ρ̃v(r). The arrow indicates the location
of rc.
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0  0.5  1  1.5  2  2.5  3
smooth
Figure 5: The AILPS Vps(r) for Si generated by the target
density ρ̃v(r) in Fig. 4.
 0  0.5  1  1.5  2  2.5  3
curate. Our proposed HQ and LQ models can be used
with any existing AILPS. However, since we want to as-
sess the performance of these models for a wide class of
atomic systems, we describe here a new method for de-
termining reasonable AILPS for general atomic systems.
These pseudopotentials will be used in all our calculations
and can be transfered to other molecular and solid state
environments, but we expect (and find in cases where
comparison can be made) little change if other reason-
able AILPS are used.
Because of the one-to-one mapping between the effec-
tive one-body potential acting on a system of N electrons
and the electron density in the ground-state configura-
tion, it is possible to obtain a unique local one-body po-
tential that generates a given target density ρ̃(r) by using
Figure 6: (Color online) Radial valence density r2ρv(r) of
the Si atom using the KS method and various models using
AILPS. Parameters used for constructing this reference sys-
tem are shown in Table IV. The arrow indicates the location
of rc. Inset: The corresponding radial total density r
2ρ(r),
which is dominated by the core component for r < rc.
 0.05
 0.15
 0.25
 0.5  1  1.5  2  2.5  3
TF1/5W
 0  0.5  1  1.5  2
a KS orbital-based method in an inverse way [23]. To con-
struct an AILPS for a given atom we separate the total
electron density ρ(r) into a “core density” ρ̃c(r), which is
supposed not to vary significantly in other molecular or
solid state environments, and the target“valence density”
ρ̃v(r) where
ρ(r) = ρ̃v(r) + ρ̃c(r). (39)
Because DFT requires only the electron density, we can
take a more general view of what is meant by the core and
valence components than is used in most orbital-based
methods. Here, we directly construct a smooth target
valence density for the Nv = N − Nc valence electrons,
with Nc chosen to be the number of electrons in the noble
gas configuration.
Our proposed target valence density ρ̃v(r) for atoms
equals the full KS density ρKS(r) outside a core of radius
rc, and is designed to be small and slowly varying inside
rc. The functional form we take is
ρ̃v(r) =
tρKS(rc) + a0r
q exp[−rp(a
r2)] r ≤ rc
ρKS(r) r > rc
Figure 4 gives an example of ρ̃v(r) for Si that will be
discussed in more detail below. We find most results are
insensitive to the details of our fitting procedure. Param-
eter values for a variety of atomic systems are given in
the Appendix in Table IV along with the physical and
technical considerations that guided our choice of this
particular form for ρ̃v(r). The Appendix also discusses
some of the general issues that arise in using these atomic
AILPS in other environments.
The local pseudopotential is directly related to the
effective one body potential that reproduces ρ̃v(r) ex-
actly when using the full KS theory. Following previous
work [24, 25], for a given ρ̃v(r), the inverse-KS equations
are solved to get the effective one-body screened poten-
tial Vscr(r). The desired ab initio local pseudopotential
Vps(r) is then obtained by subtracting the Hartree po-
tential and the exchange-correlation potential:
Vps(r) = Vscr(r)− VH(r; [ρ̃v])− Vxc(r; [ρ̃v]). (41)
This relatively expensive procedure to determine Vps(r)
requires the use of orbitals. However it needs to be done
only once for each atom, and the resulting Vps(r) can
then be used in a variety of other environments if the
atomic core densities remain essentially constant.
Once suitable Vps(r) have been determined (by this or
other means) they can be incorporated in Vext(r) in dif-
ferent ways, depending on the particular system of inter-
est. OF-DFT theory can then be used to determine the
valence density ρv(r) in direct analogy to the all-electron
calculations for the full atomic potential in Eq. (38):
▽2 + VH(r; [ρv]) + Vxc(r; [ρv])
+Vext(r) + VP (r; [ρv])}ψv(r) = µψv(r) (42)
where
ρv(r) = |ψv(r)|
2. (43)
A simple and direct test of OF-DFT is to use Eqs. (42)
and (43) for the same atomic system for which Vps(r) was
constructed. Thus we take Vext(r) = Vps(r) for a given
atom as input data in Eq. (42). The valence density
ρv(r) predicted by the OF KPs is determined from Eqs.
(42) and (43), and can be directly compared to the exact
target density ρ̃v(r) for this atomic system given by the
full KS theory.
This is illustrated in Fig. 4, which shows the input
target valence density ρ̃v(r) for Si used in the inverse-
KS process. The total KS density ρKS(r) equals ρ̃v(r)
for r ≥ rc, indicated by the arrow in Fig. 4, and then
increases rapidly for r < rc, reaching a very large value
at the nucleus, ρKS(0) = 1754. In contrast, the proposed
target valence density ρ̃v(r) remains small and relatively
slowly varying inside the core, with ρ̃v(0) = 0.005621.
Also shown in Fig. 4 are the predicted valence densities
ρv(r) for Si given by the LQ and HQ models. Because of
the relatively weak Vps(r) and slowly varying valence den-
sity, both the LQ and HQ KP models predict results very
close to those given by the exact KS treatment of the ki-
netic energy, and perform markedly better than they did
for the all-electron calculations using the full Coulomb
potential. Fig. 5 shows the corresponding AILPS gener-
ated by the inverse KS procedure. It is much smaller in
the core region than the full atomic potential and more
likely to be accurately treated by LR-based methods.
The radial valence density r2ρv(r) of the Si and Ar
atoms predicted by the various methods are shown in
Figs. 6 and 7. The consistency of our OF theory when
pseudopotentials are used is illustrated by the similar-
ity of the density predicted by the HQ and LQ models.
Figure 7: (Color online) Same as in Fig. 6 but for the Ar
atom.
 0.5  1  1.5  2  2.5  3
TF1/5W
 0  0.5  1  1.5  2
The slight deviations from the KS-DFT results for Ar in
the valence region in Fig. 7 are among the largest we en-
countered for all atoms tested, and could be due to the
relatively large number of valence electrons (8) compared
to core electrons (10). Further errors may arise from the
LFWV approximation in Eq. (33).
Once ρv(r) has been determined using OF-DFT, it can
be added to the known input core density ρ̃c(r) to obtain
the predicted total density, since the basic assumption of
our AILPS is that the core density remains unchanged
in different chemical environments. The core density is
defined in Eqs. (39) and (40). The inset in Fig. 6 shows
that the HQ and LQ treatment of the valence density
for the Si atom does not produce noticeable errors in the
total density, as expected, and the shell structure remains
in excellent agreement with the full KS-DFT calculations.
As a further test of OF-DFT we can compare the en-
ergy for the valence density given by the various methods
to the exact valence energy for the target valence density
determined by the inverse KS method. Table III gives
the valence energy values for the KS-DFT, and the LQ
and the HQ models. As can been seen, both the LQ
and HQ models give very good agreement with KS-DFT
and perform significantly better than the other models.
These results show that for this class of relatively weak
pseudopotentials the OF treatment of the KP is quite
satisfactory.
To test of the transferability of the present AILPS, we
also performed calculations for positive ions. Fig. 8 shows
that the ionization energies of various atoms calculated
using KS-DFT and the full atomic densities and those
from the KS-DFT using the valence densities with AILPS
are very similar. Therefore, the present AILPS are quite
transferable to these positive ions.
The ionization energy for models using the AILPS is
obtained by subtracting the valence energies for systems
with Nv and Nv − 1 electrons. Since the core electrons
are assumed to be unaltered in different chemical envi-
Table III: The total valence energy Ev[ρv] using the KS
method, the LQ and HQ models, the CAT model, and the
TFλW models. MAE, the mean absolute error (relative to
the KS method) of various OF models are given at the bot-
tom of their respective columns. Parameters used in Eq. (40)
for such systems are given in Table IV.
KS LQ HQ CAT TF1/5W TFW
Be −0.9914 −0.8955 −0.8950 −0.9583 −1.214 −0.7786
C −6.134 −6.080 −6.100 −6.345 −7.761 −5.266
N −11.04 −11.06 −11.09 −11.32 −13.93 −9.462
O −18.01 −18.09 −18.09 −18.19 −22.44 −15.30
Si −3.771 −3.738 −3.750 −3.869 −4.467 −3.350
P −6.474 −6.432 −6.455 −6.582 −7.385 −5.756
S −10.20 −10.10 −10.14 −10.24 −11.27 −9.023
Ar −21.37 −20.84 −20.91 −20.85 −22.40 −18.56
MAE 0.119 0.103 0.184 1.610 1.312
ronments, the atomic core energy is a constant that can-
cels here or in other similar applications to molecules and
solids. Limitations of the LQ and HQ models are more
evident here, but they do capture the overall periodicity
of the ionization energies well, and perform significantly
better than the other models.
VI. CONCLUSION
In summary, we propose two nonlocal OF KPs that
satisfy exact limits for small and large wavevector per-
turbations and reproduce the exact LR function in the
homogeneous limit. These are the same limits that sev-
eral current KEDFs are designed to satisfy. However,
because of the more local nature of the KP, it is much
easier to satisfy these conditions for the KP than for the
KEDF, and there may be other physical and technical
advantages arising from use the more local KP.
In general, there is no reason to believe that any LR-
based OF-DFT should work well for arbitrary systems
where the model potentials are far beyond the LR regime.
However, most chemical processes involve changes of va-
lence electron densities, which can often be described by
a weak AILPS. Thus the use of a LR-based OF-KP to-
gether with AILPS for atomic systems at least seems well
justified. The small and relatively slowly varying ρv(r)
also provides some justification for our use of the local
FWV kF (r) in Eq. (33).
When the AILPS is used, the valence densities given by
the LQ and HQ KPs are close to those given by the KS
method. Thus the particular integration pathway used
to get the total energy value becomes unimportant. The
simple pathway in Eq. (8) is especially useful, since no
coupling parameter integration is needed.
The proposed models are not only conceptually simple,
but also exact for a model system with a weak poten-
tial and a slowly-varying density. The appearance of the
atomic shell structure was found to be very sensitive to
the accuracy of the proposed KPs. The LR-based LQ and
Figure 8: (Color online) Ionization energies (shown in
Hartree) of the first and the second row atoms using the full
KS method, and various models using AILPS. The mean abso-
lute errors (relative to the full KS method) of various models
using AILPS are: KS (0.1 eV), LQ (1.8 eV), HQ (2.1 eV),
CAT (3.9 eV), TFW (3.1 eV), and TF1/5W (4.3 eV). Ioniza-
tion energies using the TF1/5W model are not shown in the
figure due to its relatively poor performance.
 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18
KSfull
HQ KPs give at best only qualitative indications of shell
structure for full atomic systems, though total energies
are surprisingly good. Still better results for atoms and
ions can be found by focusing on the valence density as
determined by a relatively weak AILPS. While these re-
sults seem promising, improved KPs are needed and fur-
ther investigation is required to see if these ideas can be
usefully applied to other relevant systems like molecules
and solids. Some initial results along these lines will be
reported elsewhere.
Acknowledgments
This work has been supported by the NSF Grant
CHE01-11104, and by the NSF-MRSEC at the University
of Maryland under Grant No. DMR 00-80008. J.D.C.
acknowledges the support from the UMCP Graduate
School Fellowship, the IPST Alexander Fellowship, and
the CHPH Bolck Grant Supplemental Fellowship. We
are grateful to Prof. Emily Carter and members of her
group for many helpful discussions and comments on an
earlier version of this paper.
Table IV: Parameters used in Eq. (A.1) for the target valence
density of various atoms. Here, p = q = 6, and t = 0.1
are used for all systems. The Vps(r) generated from these
parameterized ρ̃v(r) can then be used in OF-DFT.
Nv a0 a1 a2 rc
Li 1 2.983 × 10−4 0.05260 −6.560 × 10−3 2.135
Be 2 0.02655 0.7078 −0.2130 1.370
B 3 0.8032 5.855 −3.531 0.9714
C 4 11.53 30.17 −31.20 0.7429
N 5 98.80 112.8 −180.2 0.5978
O 6 592.7 338.9 −779.8 0.4981
F 7 2750 871.4 −2744 0.4256
Ne 8 1.052 × 104 1993 −8267 0.3707
Na 1 5.234 × 10−5 9.840 × 10−3 −6.610 × 10−4 2.904
Mg 2 9.805 × 10−4 0.04445 −5.002 × 10−3 2.233
Al 3 6.164 × 10−3 0.1273 −0.02056 1.861
Si 4 0.02942 0.3119 −0.06861 1.593
P 5 0.1149 0.6864 −0.1979 1.390
S 6 0.3841 1.390 −0.5100 1.231
Cl 7 1.133 2.632 −1.201 1.103
Ar 8 3.018 4.714 −2.624 0.9985
Appendix: AB INITIO LOCAL
PSEUDOPOTENTIALS
As discussed above, our proposed target valence den-
sity for atoms has the following form:
ρ̃v(r) =
tρKS(rc) + a0r
q exp[−rp(a
r2)] r ≤ rc
ρKS(r) r > rc
(A.1)
Here p and q are taken as even integers. The larger
they are, the smaller and more slowly varying is the va-
lence density near r = 0 but the sharper is the peak
near the core radius rc. As a compromise, we take here
p = q = 6, which generates relatively slowly varying local
pseudopotentials Vps(r). For applications in different en-
vironments, such as molecules or crystals, the core size rc
has to be small to maintain transferability of the atomic
core density. For this reason, we force ρ̃v(r = 0) to be
small by taking a small t. If t = 0, the strict vanish-
ing of the valence density near the nucleus would require
a very repulsive Vps, which is certainly undesirable for
the LR-based OF-DFT. However, if t is too large, there
will exist a long oscillatory tail outside the core in the
corresponding Vps(r). This is an undesirable feature for
transferability to other environments, as will be discussed
below. These two points constrain the value of t and we
use here t = 0.1 for all the atomic systems considered.
The four parameters a0, a1, a2, and rc are determined
by requiring continuity of the function ρ̃v(r) and its first
two derivatives at r = rc, and by satisfying the normal-
ization condition:
Nv = 4π
ρ̃v(r)r
2dr. (A.2)
Here we used the standard noble gas cores to determine
Figure 9: The inverse KS procedure generates very small
oscillations in the tail of rVps(r) for Si (shown here) and
other atoms. The two points where rVps(r) = −4 for Si are
r1 = 2.336 and r2 = 4.576 or r
= 2.139. The arrow indi-
cates the location of r2. To achieve good transferability (see
text), this rVps(r) is modified by setting rVps(r) = −4 when
r ≥ r2 = 4.576. See Fig. 5 for a large scale view.
-4.03
-4.02
-4.01
-3.99
-3.98
-3.97
-3.96
-3.95
 1.6  1.8  2  2.2  2.4  2.6  2.8  3  3.2  3.4
Nv = N −Nc, though other choices could in principle be
made.
To construct our local pseudopotential Vps(r) for
atoms, we first solve the KS equations for an atom with
the full Coulomb potential. With the KS density ρKS(r)
and Nc determined, this construction ensures that as
r → rc, the associated core density smoothly approaches
zero as O(| r − rc |
After generating the parameterized target valence den-
sity ρ̃v(r), the set of inverse-KS equations are solved
to obtain the corresponding one-body screened poten-
tial Vscr(r). In principle, the AILPS Vps(r) is then given
by Eq. (41). Using this, we find with an acceptable Nc
that essentially the same density profiles outside rc are
predicted by the LQ and HQ models for a wide range of
choices of p, q, and t. This would be expected if most fea-
tures of the resulting set of model potentials are within
the different regimes accurately described by the OF KP.
These results clearly show that OF-DFT can give accu-
rate results for this class of physically relevant and rela-
tively weak and slowly varying model potentials.
However, we found that the rVps(r) constructed in this
way for group III to group VIII elements deviates from
−Zv outside the core by a small and long range oscilla-
tion. For the 2s22px atoms, the maximum amplitude of
the oscillation is about 0.18. For the 3s23px atoms, it
ranges from 0.003 for Al to 0.037 for Ar. Similar oscil-
latory tails were observed by Wang and Stott [25], and
are thought to arise from the inability of a local pseu-
dopotential to represent both s and p orbitals of the cor-
responding nonlocal pseudopotentials. These oscillations
are so small in magnitude that they have almost no effect
on the density profile or energy of atoms, but they can
cause transferability problems when used for solids [53].
Following Wang and Stott’s approach, we replace Vps(r)
by −Zv/r at the largest point where rVps(r) = −Zv (see
Fig. (9)). The new Vps(r) can then be used for both
atomic and solid state calculations.
In other chemical environments, such molecules and
solids, the predetermined AILPS Vps(r) centered at each
nucleus are regarded as input data, and the valence den-
sities for other systems can be calculated by OF-DFT.
Results can be checked by using the full KS-DFT. For
example, for molecules and solids, the external potential
Vext(r) is a linear combination of the local atomic pseu-
dopotentials centered at each ion position RI :
Vext(r) =
Vps(r−RI) (A.3)
The valence density can then be determined from Eqs.
(42) and (43), and the total energy of valence electrons
Ev[ρv] = Ts[ρv] + EH [ρv] + Exc[ρv] +
ρv(r)Vext(r)dr
(A.4)
For the proposed LQ and HQ KPs, the value of Ts[ρv]
in Eq. (A.4) is determined by Eqs. (8) or (11) with ρ(r)
replaced by ρv(r).
[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864
(1964).
[2] R. G. Parr and W. Yang, Density-Functional Theory of
Atoms and Molecules, (Oxford University Press, 1989).
[3] R. M. Dreizler and E. K. U. Gross, Density Functional
Theory: An Approach to the Quantum Many Body Prob-
lem, (Springer-Verlag, Berlin, 1990).
[4] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[5] L. J. Sham and W. Kohn, Phys. Rev. 145, 561 (1966).
[6] R. A. King and N. C. Handy, Phys. Chem. Chem. Phys.
2, 5049 (2000); R. A. King and N. C. Handy, Mol. Phys.
99, 1005 (2001).
[7] J.-D. Chai and J.D. Weeks, J. Phys. Chem B 108, 6870
(2004).
[8] Y. A. Wang, N. Govind, and E. A. Carter, Phys. Rev. B
58, 13465 (1998).
[9] Y. A. Wang, N. Govind, and E. A. Carter, Phys. Rev. B
60, 16350 (1999).
[10] See e.g., Y. A. Wang and E. A. Carter, in Theoreti-
cal Methods in Condensed Phase Chemistry , edited by
S. D. Schwartz, Progress in Theoretical Chemistry and
Physics, (Kluwer, Boston, 2000) p. 117, and references
therein.
[11] B. Zhou, V. L. Ligneres, and E. A. Carter, J. Chem.
Phys. 122, 044103 (2005).
[12] E. Chacón, J. E. Alvarellos, and P. Tarazona, Phys. Rev.
B 32, 7868 (1985).
[13] P. Garćıa-González, J. E. Alvarellos, and E. Chacón,
Phys. Rev. A 54, 1897 (1996).
[14] P. Garćıa-González, J. E. Alvarellos, and E. Chacón,
Phys. Rev. A 57, 4857 (1998).
[15] M. Pearson, E. Smargiassi, and P. A. Madden, J. Physics:
Condens. Matter 5, 3221 (1993); E. Smargiassi and P. A.
Madden, Phys. Rev. B 49, 5220 (1994); M. Foley and P.
A. Madden, Phys. Rev. B 53, 10589 (1996).
[16] F. Perrot, J. Phys.: Condens. Matter 6, 431 (1994).
[17] J. A. Alonso and L. A. Girifalco, Phys. Rev. B 17, 3735
(1978); M. D. Glossman, L. C. Balbás, and J. A. Alonso,
Chem. Phys. 196, 455 (1995).
[18] C. Herring, Phys. Rev. A 34, 2614 (1986).
[19] L.-W. Wang and M. P. Teter, Phys. Rev. B 45, 13196
(1992).
[20] N. Choly and E. Kaxiras, Solid State Commun. 121, 281
(2002).
[21] L. H. Thomas, Proc. Cambridge Phil. Soc. 23, 542
(1927); E. Fermi, Z. Phys. 48, 73 (1928).
[22] See e.g., R. Evans, in Fundamentals of Inhomogeneous
Fluids, edited by D. Henderson (Dekker, New York,
1992), p. 85.
[23] Y. Wang and R. G. Parr, Phys. Rev.A 47, R1591 (1993).
[24] B. Zhou, Y.A. Wang, and E.A. Carter, Phys. Rev. B 69,
125109 (2004).
[25] B. Wang and M. J. Stott, Phys. Rev. B 68, 195102
(2003).
[26] A. Holas and N. H. March, Phys. Rev. A 66, 066501
(2002); I. Lindgren and S. Salomonson, Phys. Rev. A
67, 56501 (2003).
[27] L. R. Pratt, G. G. Hoffman, and R. A. Harris, J. Chem.
Phys. 88, 1818 (1988); L. R. Pratt, G. G. Hoffman, and
R. A. Harris, J. Chem. Phys. 92, 6687 (1990); G. G.
Hoffman and L. R. Pratt, Mol. Phys. 82, 245 (1994).
[28] Y.-G. Chen and J. D. Weeks, J. Chem. Phys. 118, 7944
(2003)
[29] J. P. Perdew and S. Kurth, in Density Functionals: The-
ory and Applications, edited by D. Joubert, Lecture
Notes in Physics, (Springer-Verlag, Berlin, 1998) p. 8,
and references therein.
[30] L. J. Sham, Phys Rev. A 1, 969 (1970).
[31] M. Levy and J. P. Perdew, Phys. Rev.A 32, 2010 (1985).
[32] F. W. Averill and G. S. Painter, Phys. Rev. B 24, 6795
(1981).
[33] W. Yang, Phys Rev. A 34, 4575 (1986).
[34] C. F. von Weizsäcker, Z. Physik 96, 431 (1935).
[35] J. Lindhard, K. Dan. Vidensk. Selsk. Mat. Fys. Medd.
28, 8 (1954).
[36] W. Jones and W. H. Young, J. Phys. C 4, 1322 (1971).
[37] Y. Tomishima and K. Yonei, J. Phys. Soc. Jpn. 21, 142
(1966).
[38] E. H. Lieb, Rev. Mod. Phys. 53, 603 (1981).
[39] G. K.-L. Chan, A. J. Cohen, and N. C. Handy, J. Chem.
Phys. 114, 631 (2001).
[40] P. K. Acharya, L. J. Bartolotti, S. B. Sears, and R. G.
Parr, Proc. Natl. Acad. Sci. USA 77, 6978 (1980).
[41] J. L. Gázquez and J. Robles, J. Chem. Phys. 76, 1467
(1982).
[42] L. J. Bartolotti and P. K. Acharya, J. Chem. Phys. 77,
4576 (1982).
[43] P. K. Acharya, J. Chem. Phys. 78, 2101 (1983).
[44] D. J. Singh, Phys. Rev. B 48, 14099 (1993).
[45] J.-D. Chai, Ph.D. dissertation, University of Maryland,
2005.
[46] Instead of solving the second-order differential equation
numerically for the corresponding weight function ω(x)
(see Eq. (14) of Ref. [13]), we used the parametrized form
in the appendix of Ref. [13]. Although the parametriza-
tion was made for β = −1/2 (note the sign change of β
in Ref. [13]), the β dependence of ω(x) is small for small
β, as argued in Ref. [14]. Our results for the CAT model
using this parametrized ω(x) are very close to their re-
ported results.
[47] P. A. M. Dirac, Proc. Cambridge Phil. Soc. 26, 376
(1930).
[48] D. M. Ceperley, Phys. Rev. B 18, 3126 (1978); D. M.
Ceperley and B.J. Alder, Phys. Rev. Lett. 45 , 566
(1980).
[49] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048
(1981).
[50] N. H. March, Phys. Lett. A 113, 476 (1986).
[51] A. Holas and N. H. March, Phys. Rev.A 44, 5521 (1991).
[52] A. M. Abrahams and S. L. Shapiro, Phys. Rev. A 42,
2530 (1990).
[53] D. M. Bylander and L. Kleinman, Phys. Rev. B 55, 9432
(1997).

  In the density functional (DF) theory of Kohn and Sham, the kinetic energy of
the ground state of a system of noninteracting electrons in a general external
field is calculated using a set of orbitals. Orbital free methods attempt to
calculate this directly from the electron density by approximating the
universal but unknown kinetic energy density functional. However simple local
approximations are inaccurate and it has proved very difficult to devise
generally accurate nonlocal approximations. We focus instead on the kinetic
potential, the functional derivative of the kinetic energy DF, which appears in
the Euler equation for the electron density. We argue that the kinetic
potential is more local and more amenable to simple physically motivated
approximations in many relevant cases, and describe two pathways by which the
value of the kinetic energy can be efficiently calculated. We propose two
nonlocal orbital free kinetic potentials that reduce to known exact forms for
both slowly varying and rapidly varying perturbations and also reproduce exact
results for the linear response of the density of the homogeneous system to
small perturbations. A simple and systematic approach for generating accurate
and weak ab-initio local pseudopotentials which produce a smooth slowly varying
valence component of the electron density is proposed for use in orbital free
DF calculations of molecules and solids. The use of these local
pseudopotentials further minimizes the possible errors from the kinetic
potentials. Our theory yields results for the total energies and ionization
energies of atoms, and for the shell structure in the atomic radial density
profiles that are in very good agreement with calculations using the full
Kohn-Sham theory.

Orbital-Free Density Functional Theory: Kinetic Potentials and Ab Initio Local
Pseudopotentials
Jeng-Da Chai∗1 and John D. Weeks1, 2
Institute for Physical Science and Technology,
and Department of Chemistry and Biochemistry,
University of Maryland, College Park, MD 20742
(Dated: November 1, 2018)
In the density functional (DF) theory of Kohn and Sham, the kinetic energy of the ground state
of a system of noninteracting electrons in a general external field is calculated using a set of orbitals.
Orbital free methods attempt to calculate this directly from the electron density by approximating
the universal but unknown kinetic energy density functional. However simple local approximations
are inaccurate and it has proved very difficult to devise generally accurate nonlocal approximations.
We focus instead on the kinetic potential, the functional derivative of the kinetic energy DF, which
appears in the Euler equation for the electron density. We argue that the kinetic potential is more
local and more amenable to simple physically motivated approximations in many relevant cases, and
describe two pathways by which the value of the kinetic energy can be efficiently calculated. We
propose two nonlocal orbital free kinetic potentials that reduce to known exact forms for both slowly
varying and rapidly varying perturbations and also reproduce exact results for the linear response of
the density of the homogeneous system to small perturbations. A simple and systematic approach
for generating accurate and weak ab initio local pseudopotentials which produce a smooth slowly
varying valence component of the electron density is proposed for use in orbital free DF calculations
of molecules and solids. The use of these local pseudopotentials further minimizes the possible errors
from the kinetic potentials. Our theory yields results for the total energies and ionization energies
of atoms, and for the shell structure in the atomic radial density profiles that are in very good
agreement with calculations using the full Kohn-Sham theory.
I. INTRODUCTION
Density-functional theory (DFT) has become one of
the most powerful tools for investigating the electronic
structure of large complex systems. In principle, as
shown by Hohenberg and Kohn [1], the exact ground
state energy of a system of N electrons can be formally
written as a functional E[ρ] of only the electron density
ρ(r), a function of three variables, and the external field
Vext(r). Determining the energy and other ground state
properties from such an approach could dramatically
reduce the computational cost for large systems when
compared with traditional quantum chemistry methods,
which deal with wavefunctions involving coordinates of
all N electrons [2, 3].
Kohn and Sham (KS) [4, 5] showed that E[ρ] can be
usefully partitioned into the following set of terms:
E[ρ] = Ts[ρ] + EH [ρ] + Exc[ρ] +
ρ(r)Vext(r)dr. (1)
Here Ts[ρ] is the noninteracting kinetic energy density
functional (KEDF), which gives the kinetic energy of
a model system of N noninteracting electrons in a self-
consistent field chosen so that the ground state density
∗E-mail: jdchai@berkeley.edu. Present address: Molecular
Foundry, Materials Sciences Division, Lawrence Berkeley National
Laboratory, and Department of Chemistry, University of Califor-
nia, Berkeley, California 94720
equals ρ(r),
EH [ρ] ≡
ρ(r)ρ(r′)
|r− r′|
drdr′ (2)
is the classical electron-electron potential energy (Hartree
energy) and Exc[ρ] is the exchange-correlation energy (in-
cluding the difference between the interacting and non-
interacting kinetic energy and the difference between
the quantum and classical electron-electron potential en-
ergy). The last term on the right of Eq. (1) is the only
term that depends explicitly on the external potential
Vext(r). Atomic units are used throughout the paper.
If all these functionals were known, then the density
ρ(r) could be obtained from the variational principle (Eu-
ler equation) associated with minimizing Eq. (1):
µ = VTs(r; [ρ]) + Veff (r; [ρ]), (3)
and the total energy of the inhomogeneous system could
then be determined from the energy functional E[ρ]. All
other physical quantities related to the ground-state den-
sity could also be computed. Here µ is the chemical po-
tential (the Lagrange multiplier associated with the nor-
malization condition
ρ(r)dr = N), and Veff (r; [ρ]) is
an effective one-body potential defined by
Veff (r; [ρ]) ≡
δρ(r)
EH [ρ] + Exc[ρ] +
ρ(r)Vext(r)dr
= VH(r; [ρ]) + Vxc(r; [ρ]) + Vext(r), (4)
where
VH(r; [ρ]) ≡ δEH [ρ]/δρ(r) =
ρ(r′)
|r− r′|
dr′ (5)
http://arxiv.org/abs/0704.1878v1
is the Hartree potential, and Vxc(r; [ρ]) ≡ δExc[ρ]/δρ(r) is
the exchange-correlation potential. Similarly we interpret
VTs(r; [ρ]) ≡ δTs[ρ]/δρ(r) (6)
as the kinetic potential (KP) arising from the KEDF [6].
Further progress requires an accurate determination
of the noninteracting kinetic energy, whose magnitude is
much larger than the exchange-correlation energy. The
initial development of DFT as a practical computational
method was made possible by KS’s realization that the
numerical value of the noninteracting kinetic energy can
be exactly calculated, not directly from the density itself
using Ts[ρ], but by introducing a set of N one-electron
wave functions (orbitals) satisfying the N coupled KS
equations that describe the model system [4, 5].
Research could then focus on determining the remain-
ing small term Exc[ρ]. Here even local density approx-
imations have often proved useful. Through the efforts
of many workers we now have generally accurate expres-
sions for Exc[ρ]. Using these along with the KS orbitals
to calculate the kinetic energy, one can accurately cal-
culate both the total energy E[ρ] and the ground-state
density ρ(r) for a wide variety of systems.
However, the use of the KS orbitals usually generates
a relatively expensive O(N3) scaling of computational
cost with the number of electrons. While this scaling is
much better than that of most standard methods that
include correlation energy, calculations for large systems
remain problematic. This remaining bottleneck could be
removed if there were an accurate treatment of the kinetic
energy in terms of the electron density only [2, 3, 7].
To that end there has been considerable effort in-
vested in developing “orbital-free”density functional the-
ory (OF-DFT) by making direct approximations for Ts[ρ]
[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. While
earlier simple local density approximations for Ts[ρ] like
those used in the Thomas-Fermi (TF) model [21] are very
inaccurate, there have been two main advances in recent
work that offer prospects for significant improvements.
The first is the introduction of nonlocal KEDFs that
reproduce known exact results for very slowly varying or
very rapidly varying fields and give the exact linear re-
sponse (LR) of the density of the uniform model system
to small perturbations. Similar ideas have been success-
fully applied to classical nonuniform fluids [22]. The sec-
ond advance is to focus not on the total density but on
the smaller and more slowly varying valence electron den-
sity as described by a weak pseudopotential acting only
on the valence electrons. While conventional pseudopo-
tential methods use orbitals, recently developed ab ini-
tio local pseudopotential (AILPS) methods determine the
unique local one-body potential producing a given target
valence density by solving the KS equations inversely, us-
ing the one-to-one mapping between density and poten-
tials in DFT [23]. For OF-DFT with LR-based KEDF’s,
the use of pseudopotentials not only can reduce the com-
putational cost, but also can improve its accuracy, since
the system will be closer to the LR regime where Ts[ρ] is
designed to be accurate. [24, 25]. Indeed, very promising
results using such OF-DFT methods have been obtained
for a variety of nearly free-electron-like metals.
However, existing KEDF’s have not yet achieved chem-
ical accuracy for systems with localized and more rapidly
varying electron densities like molecules or for covalent
or ionic solids. The main problem is that the exact Ts[ρ]
is highly nonlocal, and we have little idea of the func-
tional form of the nonlocality for densities far from the
LR regime. It has proved very difficult to understand
what errors an approximate nonlocal Ts[ρ] will produce
in the density as determined by the Euler equation with
a general Vext(r).
We explore here a different way to attack this basic
problem. The exact Ts[ρ] can be formally obtained from
VTs(r; [ρ]) by functional integration over density changes
in all regions of space [6, 18]. Because of this integration
Ts[ρ] is a more nonlocal functional of the density than
is VTs(r; [ρ]). More detailed arguments arriving at this
same conclusion have been recently presented [26]. Since
most problems in devising accurate approximations for
Ts[ρ] have arisen from the nonlocality, this suggests it
could be worthwhile to try to develop approximations
for the KP VTs(r; [ρ]) itself.
To illustrate this point, Chai and Weeks [7] added a
simple gradient correction to the original local TF KP
for atoms [21], with a coefficient chosen to reproduce
the exact boundary condition of exponential decay of the
electron density far from the nucleus. Though quanti-
tative results were not obtained, the resulting modified
Thomas-Fermi (MTF) model gave energies for atoms and
for closed-shell diatomic molecules that showed notable
improvements when compared to the original TF and re-
lated gradient corrected KEDF models. However, the
local gradient correction used in the MTF KP cannot re-
produce the oscillatory atomic shell structure and it does
not satisfy the exact LR behavior in the homogeneous
limit. It is clear that nonlocality even in the KP must be
taken into account to achieve more accurate results.
We propose here new nonlocal approximations for the
KP using ideas similar to those employed for the nonlocal
KEDFs. These new KPs satisfy the exact LR condition
in the uniform limit, and reproduce known exact limiting
forms of VTs(r; [ρ]) both for very slowly varying and very
rapidly varying perturbations. As will become clear, the
nonlocality in our KP is determined by the requirement
that LR is exactly satisfied, and it is much easier to en-
sure that LR holds for the KP than it is for analogous
KEDF models. We believe this level of nonlocality in the
KP may suffice in many cases when used in conjunction
with AILPS methods to describe slowly varying valence
density components closer to the LR regime.
The remainder of this paper is organized as follow.
Section II will discuss some general pathways connect-
ing Ts[ρ] and VTs(r; [ρ]). Section III will describe limit-
ing forms of the KEDF and KP for slowly varying and
rapidly varying perturbations, and discuss LR theory, an
exact theory for the response of density of the uniform
electron gas to small perturbations. Section IV will de-
velop two nonlocal KPs incorporating both the correct
limiting forms of the exact KP and the exact LR of the
free-electron gas. Section V will compare the numerical
results of the present method for atoms with the KS-DFT
and other KEDFs, both for all-electron calculations and
for valence electrons using the AILPS. We find that the
use of AILPS indeed reduces errors arising from nonlocal-
ity in these approximate KPs or KEDFs, which give very
accurate results for the relatively slowly varying valence
densities. Our conclusions are given in Section VI.
II. PATHWAYS FROM VTs(r; [ρ]) TO Ts[ρ]
If Ts[ρ] is known, VTs(r; [ρ]) can be simply computed
by functional differentiation. However, there is no unique
way of determining Ts[ρ] from a given VTs(r; [ρ]). Many
possible pathways can be used to construct Ts[ρ] by func-
tional integration of VTs(r; [ρ]) [27, 28]. If the exact
VTs(r; [ρ]) is used and the integration is carried out ex-
actly, then all pathways would give the same exact result
for Ts[ρ]. However, when an approximate VTs(r; [ρ]) is
used, different pathways will give different results for the
kinetic energy. But this “thermodynamic inconsistency”
is small if reasonably good approximations are used, since
the integration tends to smooth out local errors that may
exist in the density [7, 28].
More problematic is the fact that most pathways re-
quire additional results for partially coupled systems as
the external field or density perturbation is gradually
turned on, which adds to the computational burden. In
particular, most earlier work has used a “potential en-
ergy pathway”, where the external potential is scaled by
a coupling parameter [7, 27]. The kinetic energy can then
be found by subtracting the potential energy (calculated
from the potential energy density functionals) from the
total energy. However, this pathway is expensive, since
one has to solve the Euler equation (3) for each partially
coupled V λext(r) (with the same µ), to determine the cor-
responding ρλ.
A. Herring’s Pathway
However, Herring showed there is a particular path-
way arising from exact scaling relations between the non-
interacting kinetic energy Ts[ρ] with respect to the coor-
dinate r in ρ(r) where very simple results involving only
the final density can sometimes be found [6, 18]. If the
coordinate r is scaled to αr, the normalized scaled den-
sity is ρα(r) = α
3ρ(αr). It is easy to show that the exact
Ts[ρ] then obeys
α2Ts[ρ] = Ts[ρα]. (7)
For isolated systems, such as atoms and molecules, the
density and its derivatives to all order vanish far from the
nuclei. For such systems, when Eq. (7) is differentiated
with respect to α, and the partial derivative is evaluated
at α = 1, we find the formally exact result
Ts[ρ] =
VTs(r; [ρ])∇ · (rρ(r))dr. (8)
Therefore, once the kinetic potential VTs(r; [ρ]) is
known for some given ρ(r), the numerical value of Ts[ρ]
can then be immediately determined from Eq. (8). Since
there is no need to perform a coupling parameter integra-
tion over the change of density or potential, this scheme
is not only fast, but also numerically reliable. The final
form of Eq. (8) is essentially the viral theorem, and is
directly related to the force on molecules [29, 30, 31, 32].
Note that this simple and exact pathway holds only
for the noninteracting Ts[ρ] [31], which again shows the
virtues of the KS partitioning of the total energy. We will
use Eq. (8) as the basic pathway to determine the numer-
ical value of Ts[ρ] from a given approximate VTs(r; [ρ])
for most calculations in this paper. However, Eq. (8)
does not hold for extended solid-state systems because of
nonvanishing boundary terms, and thus far we have not
found an exact and simple way of including them.
Fortunately, there is another class of computationally
efficient“density pathways”that can be used for extended
systems, as we now show. Density pathways can also
be used for atomic and molecular systems to check the
accuracy of the VTs used, since results using the exact
VTs would be independent of path [27, 28].
B. Density Pathways
The change in the kinetic energy can be formally re-
lated to a coupling parameter integration, where the den-
sity changes from some known value at λ = 0 to the final
density at λ = 1:
Ts[ρ] = Tλ=0 +
drVTs(r; [ρλ])
∂ρλ(r)
In most cases a simple linear density pathway will suffice.
Here the density ρ(r) is linearly scaled by a coupling pa-
rameter λ from some uniform reference density ρ0 natu-
rally chosen to be the uniform electron density N/V in
extended systems:
ρλ(r) = ρ0 + λ[ρ(r) − ρ0]. (10)
Then Eq. (9) becomes
Ts[ρ] = Tλ=0 +
drVTs(r; [ρλ])[ρ(r) − ρ0] (11)
Here Tλ=0 is the kinetic energy of the uniform system,
i.e., the Thomas-Fermi kinetic energy TTF [ρ0]. For ex-
tended systems, where the Herring’s pathway cannot be
used, this density pathway appears to be a good way to
compute T . Other density pathways, like the square-root
pathway introduced by Chen and Weeks [28] to describe
nonuniform hard sphere fluids, can be defined, and have
proved useful in certain applications, but we do not con-
sider them here.
Note from Eq. (10) that ρλ(r) depends only on the fi-
nal density, so evaluation of Eq. (11) is straightforward
and this pathway is computationally efficient. Unlike the
potential energy pathway where the external potential is
scaled, there is no need to solve the Euler equation (3)
for its corresponding external potential V λext(r) at each λ.
However, for isolated systems, where ρ0 = 0, this path-
way is likely to be less accurate than Herring’s pathway,
since it does not automatically satisfy the virial theorem.
III. EXACT LIMITS AND LINEAR RESPONSE
THEORY
Although the exact Ts[ρ] is still unknown, several limit-
ing forms have been discovered for particular density dis-
tributions. These provide important cornerstones that
can be used to construct accurate KEDFs and KPs in
many cases, as will be seen below.
In particular, the Thomas-Fermi (TF) KEDF [21] is
known to be exact for a uniform system:
TTF [ρ] = CF
ρ5/3(r)dr, (12)
where CF =
(3π2)2/3. The TF KEDF TTF [ρ] is de-
rived by local use of uniform free-electron gas model, and
is exact for a system with an infinite number of electrons.
The corresponding expression for the TF KP is
VTF (r; [ρ]) ≡ δTTF [ρ]/δρ(r) =
2/3(r). (13)
This depends only on the local value of ρ2/3(r) and thus
formally is more local than the TF KEDF, whose func-
tional dependence on ρ involves the density at all r. Of
course in this simple case the functional integration of
Eq. (13) can be carried out exactly to yield Eq. (12),
but this cannot be done in general and the nonlocality of
Ts[ρ] has proved problematic.
Results for nonuniform systems are best described in
Fourier space. For a very slowly varying perturbation of
the density, the second-order gradient expansion is exact
[33]. It is easy to see that results correct to second order
at small wavevectors are given by
TTF 1
W [ρ] ≡ TTF [ρ] +
TW [ρ], (14)
where
TW [ρ] ≡
|▽ρ(r)|2
dr (15)
is the von Weizsäcker (W) KEDF [34].
TW [ρ] is exact for a system with one or two electrons,
or where the density can be accurately described by a sin-
gle orbital. Moreover it has been argued [8, 9, 10, 12] that
TW [ρ] gives the correct leading order term for a rapidly
varying perturbation with only high wavevector compo-
nents and that the next order correction is reproduced
TW− 3
TF [ρ] ≡ TW [ρ]−
TTF [ρ]. (16)
The W KP is
VW (r; [ρ]) ≡ δTW [ρ]/δρ(r) =
|▽ρ(r)|2
ρ2(r)
▽2ρ(r)
If we represent the full density by an effective single or-
bital function ψ(r),
ρ(r) = |ψ(r)|2 (18)
then the W KP can be written in a compact form that
will later prove useful:
VW (r; [ρ]) = −
▽2ψ(r)
2ψ(r)
. (19)
Finally, the linear response of the density of a uniform
non-interacting electron gas with density ρ0 to a small
perturbation δV (k) = ǫke
ik·r is exactly known [35],
δρ(k) = χL(q)δV (k). (20)
q ≡ k/2kF (21)
is a dimensionless wavevector, where
kF ≡ (3π
1/3 (22)
is the Fermi wavevector (FWV) and k ≡ |k|. The LR
function χL(q) has the form
χL(q) = −
F−1L (q)
1− q2
1 + q
, (23)
where
FL(q) ≡
1− q2
1 + q
has been called the Lindhard function [10].
It is known that the weak logarithmic singularity at
q = 1 in F−1L (q) is responsible for Friedel oscillations, and
may also be important for the appearance of atomic shell
structure. This singularity further divides the Lindhard
function into two branches in Fourier space: the low-
momentum (q < 1 ) or the low-q (LQ) branch, and the
high-momentum (q > 1) or the high-q (HQ) branch [10].
Figure 1: Linear response functions of a uniform system of
noninteracting Fermions as given by the TF, W, and MTF
(see Ref. [7]) models.
 0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2  2.2  2.4  2.6  2.8  3
exact
The dimensionless response function arising from the
TF KEDF is FTF (q) = 1, and that from the W KEDF is
FW (q) = 3q
2 [36]. Clearly, no linear combination of the
TF and the W KEDFs can reproduce the exact Lindhard
function in Eq. (24). This has the following two limits
[10],
FL(q) =
1 + q2/3 +O(q4) q ≪ 1
3q2 − 3/5 +O(q−2) q ≫ 1
It should be noted that the expansions for both the
low-q and high-q limits are correct to all orders in per-
turbation theory, but valid only in the appropriate limits
in Fourier space. On the other hand, the LR theory is
valid for all wavevectors, but is only accurate for small
perturbations. Therefore, the regime where the response
functions of the two limiting KEDFs deviate from the
exact LR function gives an indication of the range of
wavevectors where the two limiting forms are inaccurate.
As shown in Fig. 1, the response function F−1TF (q) has
no momentum dependence and is only exact at q = 0.
The response function F−1W (q) is exact asymptotically at
high q, and remains fairly accurate for q & 2, but is di-
vergent in the low-q branch, and fails completely for the
nearly uniform electron gas. In contrast, the MTF model
[7] gives a reasonably accurate average description of the
exact response function, especially in the important re-
gion near the singularity at q = 1.
IV. CONSTRUCTION OF NONLOCAL
KINETIC POTENTIALS
A. Kinetic Energy Density Functionals TTFλW [ρ]
and TWλTF [ρ]
Simple linear combinations of the two limiting KEDF’s
in Eqs. (12) and (15), such as the TFλW KEDF [33, 37,
38, 39]
TTFλW [ρ] ≡ TTF [ρ] + λTw[ρ] (26)
and the WλTF KEDF [40, 41, 42, 43]
TWλTF [ρ] ≡ TW [ρ] + λTTF [ρ] (27)
have been widely studied for several decades. The value
of the parameter λ was either determined empirically for
getting good atomic energy or obtained by some semi-
classical arguments.
The advantage of these approaches is the ability to gen-
erate a family of simple KEDF’s easily. It has been shown
empirically that the TF1/5W model can give good val-
ues for atomic energies, but the predicted density profiles
are generally not very accurate, both near and far away
from the nucleus. The TTFλW [ρ] and TWλTF [ρ] function-
als give the correct leading term in the density response
to a slowly-varying perturbation and a rapidly-varying
perturbation respectively, and with particular choices of
λ as in Eqs. (14) and (16) they can reproduce the next
order term. Unfortunately, they then will have an in-
correct leading term in the opposite limit, unless λ = 1.
However it has been shown that TTFW [ρ] with λ = 1
always overestimates the exact Ts[ρ] for various systems
[10]. Finally, none of these functionals can reproduce the
exact response function FL(q) in the homogeneous limit.
Since these models fail to satisfy all the known limits,
and nonlocality in Ts[ρ] is not correctly described, it is
also not surprising that atomic shell structure is missing
in these approaches.
B. Combining TF and W Kinetic Potentials
We argue that it may be more profitable to take advan-
tage of known limiting forms of the KP, rather than the
KEDF, and develop approximations for the more local
VTs(r; [ρ]) directly. Again we can rely on known results
in the linear response regime when the density variations
are not too large.
From Eqs. (14) and (16), the following linear combina-
tions of the TF KP and the W KP in Eqs. (13) and (17)
can reproduce exact results to second order for very small
and very large wavevector perturbations respectively:
VTs(r; [ρ]) ≈
VTF (r; [ρ]) +
VW (r; [ρ]) q ≪ 1
VW (r; [ρ])−
VTF (r; [ρ]) q ≫ 1
Since VTF (r; [ρ]) and VW (r; [ρ]) are the only components
up to second order of the two exact limiting forms of the
KP, we can combine them in analogy to the TFλW and
WλTF models and arrive at generalized KPs.
However, instead of combining them using a fixed pa-
rameter λ, it seems natural to represent them in Fourier
space and allow a wavevector dependence in λ = λ(q) to
connect the limiting forms. The λ(q) can then be chosen
in a very simple way so that the exact LR function is
Figure 2: Weight function f̂(q) for the HQ and LQ KPs.
 0  0.5  1  1.5  2  2.5  3
reproduced for a uniform system with density ρ0. In this
way the LR function bridges the exact limits at large and
small wavevectors, and if the theory is applied to weak
perturbations in the linear response regime for interme-
diate wavevectors we can expect very accurate results.
Here, we derive such generalized KPs based on the KP
for the WλTF model.
C. HQ Kinetic Potential
In analogy to the WλTF model in Eq. (27) we look for
a kinetic potential of the form
V 0HQ(k) = VW (k) + λHQ(q) VTF (k)
= VTF (k) + VW (k) + f̂(q) VTF (k), (29)
where
q = k/2kF (30)
is a dimensionless wavevector normalized by the FWV kF
in Eq. (22) of a uniform reference system with density
ρ0 and f̂(q) = λHQ(q) − 1. The superscript 0 in V
indicates use of a uniform reference system. For a small
perturbation, we can linearize the V 0HQ(k) in Eq. (29).
Requiring that it satisfy LR exactly then determines the
weight function f̂(q) as
f̂(q) = FL(q)− 3q
2 − 1. (31)
See Fig. 2.
We refer to Eq. (29) with Eq. (31) as the HQ KPmodel.
It reproduces the correct high-q limit in Eq. (28) up to the
second order. However, unlike Eq. (16), it also satisfies
the correct low-q limit to leading order and gives exact
results for all q in the linear response regime. Inverse
Fourier transform of Eq. (29) then gives
V 0HQ(r; [ρ], kF ) = VTF (r; [ρ]) + VW (r; [ρ])
f(|r− r′|; kF ) VTF (r
′; [ρ])dr′.
This expression is directly useful for extended systems
where a reasonable ρ0 can be defined. For isolated sys-
tems such as atoms and molecules where the density van-
ishes far from the nuclei, it seems natural to replace kF
in Eq. (32) by the local Fermi wavevector (LFWV)
kF (r) ≡ (3π
2ρ(r))1/3, (33)
though errors may be introduced for rapidly varying den-
sity distributions. Using Eq. (13), this yields the general
form of our proposed HQ kinetic potential:
VHQ(r; [ρ], kF (r)) = VTF (r; [ρ]) + VW (r; [ρ])
f(|r− r′|; kF (r)) ρ
2/3(r′)dr′.
Note that the last term in Eq. (34) is most easily com-
puted in Fourier space as
3(2π)3
f̂(k/2kF (r))ρ
2/3(k)e−ik·rdk, (35)
As can been seen in Eq. (31), the weight function
f̂(k/2kF (r)) is determined analytically. Unlike the LR-
based KEDF approaches, no first-order differential equa-
tion is needed to solve for the weight function in Fourier
space. For atomic systems where the LFWV is used, the
convolution in Eq. (35) must be carried out numerically,
which will lead to a quadratic scaling of the HQ model
(and related LQ model described below) in the numberN
of electrons. For extended systems where one can expand
about the local density ρ0, one can use fast Fourier trans-
forms (FFT’s) for a much more efficient computation of
this integral [44].
D. LQ Kinetic Potential
In analogy to the TFλW model in Eq. (26), we could
similarly generate a KP that is accurate to second order
at low q while still reproducing the leading term at high
q. However this is numerically less useful because the
analogue of Eq. (35) involves the Fourier transform of
VW (r
′; [ρ]), which cannot be simply expressed in terms
of the density. Instead, by empirically taking a properly
chosen component of the density outside the integral we
find that Eq. (34) can be modified to produce a new LQ
KP that is accurate to second order at low q and first
order at high q [45]:
VLQ(r; [ρ], kF (r)) = VTF (r; [ρ]) + VW (r; [ρ])
1/6(r)
f(|r− r′|; kF (r)) ρ
1/2(r′)dr′
Extending these ideas we have constructed a modified
KP that satisfies LR everywhere and is accurate to sec-
ond order at both low and high q [45]. However the
functional form is much more complicated, and little ad-
ditional accuracy is gained from the improved behavior
at very small or very large wavevectors, since all forms
use LR to interpolate for intermediate wavevectors, and
this is where most errors arise in practice. Thus we will
report results here only for the HQ and LQ models.
V. RESULTS FOR ATOMS USING THE HQ
AND LQ MODELS
For completeness and to compare to earlier work, we
first briefly discuss all-electron calculations using the pro-
posed HQ and LQ model KPs and the full atomic poten-
tials. We then describe results using AILPS methods.
These are compared with the KS-DFT, the TFλW mod-
els, and the CAT model introduced by Chacón, Alvarel-
los, and Tarazona [12, 13, 14]. The CAT model is a
LR-based KEDF method, which gives some indication of
shell structure. We employed the latest version, which
uses a nonlocal two-body Fermi wave vector with a pre-
scribed functional form depending on an empirical pa-
rameter β = 1/2 (defined in Eq. (3) of Ref. [14]). This
caused the numerical calculations [46] to be consider-
ably more costly than those of the LQ or HQ models,
which used the local Fermi wave vector as in Eq. (33).
All calculations are spin-restricted and use the local den-
sity approximation (LDA) [47, 48, 49] for the exchange-
correlation functional.
A. All-electron calculations
All-electron calculations consider the density response
to the large and rapidly varying nuclear potential. Since
the system is far from the linear response regime, quanti-
tative results from the HQ and LQ models (or from LR-
based KEDF methods) cannot be expected. However, by
incorporating exact results for very large and very small
wavevectors, these models do correct major deficiencies
of the purely local TF model (which, e.g., predicts an
infinite density at the nucleus!) and even give some qual-
itative indications of atomic shell structure.
The numerical method use the Pauli kinetic potential
V P (r; [ρ]) [50, 51], defined as
V P (r; [ρ]) ≡ VTs(r; [ρ]) − VW (r; [ρ]). (37)
Since VW (r; [ρ]) is the exact KP for a system where the
density can be accurately described by a single orbital,
if V P (r; [ρ]) is omitted, one would essentially obtain the
ground state density of the corresponding Boson system,
where all the electrons are in the same orbital. If we
represent the full density by a single orbital function ψ(r),
so that ρ(r) and VW (r; [ρ]) can be written in the forms of
Eq. (18) and Eq. (19) respectively, we can then combine
V P (r; [ρ]) with the one-body potential Veff (r; [ρ]) in Eq.
(3), and derive a Schrödinger-like equation for the Bose
orbital ψ(r),
▽2 + Veff (r; [ρ]) + V
P (r; [ρ])
ψ(r) = µψ(r). (38)
Table I: Atomic energy E using the KS, LQ, HQ, CAT and
the TFλW models in all-electron calculations. MAPE, the
mean absolute percentage error (relative to the KS method) of
various OF models are given at the bottom of their respective
columns.
KS LQ HQ CAT TF1/5W TFW
He −2.834 −2.565 −2.437 −2.675 −2.911 −1.559
Ne −128.2 −134.3 −126.6 −126.2 −129.5 −86.40
Ar −525.9 −545.9 −512.2 −515.1 −526.2 −375.5
Kr −2750 −2805 −2621 −2712 −2748 −2099
Xe −7229 −7306 −6844 −7141 −7214 −5701
Be −14.45 −14.39 −13.64 −14.11 −14.71 −8.699
Mg −199.1 −207.9 −195.7 −195.2 −200.0 −136.4
C −37.42 −38.97 −36.85 −37.25 −38.41 −24.01
N −54.02 −56.71 −53.59 −53.84 −55.39 −35.33
O −74.47 −78.39 −74.02 −74.08 −76.11 −49.39
Si −288.2 −300.4 −282.5 −282.2 −288.9 −200.5
P −339.9 −354.1 −332.7 −332.8 −340.6 −238.3
S −396.7 −412.8 −387.7 −388.5 −397.3 −279.9
MAPE 4.06% 3.42% 1.82% 1.10% 32.0%
Table II: Electron density at the nucleus ρ(0), using the KS,
LQ, HQ, CAT and the TFλW models in all-electron calcula-
tions. MAPE, the mean absolute percentage error (relative to
the KS method) of various OF models are given at the bottom
of their respective columns.
KS LQ HQ CAT TF1/5W TFW
He 3.525 3.088 2.742 3.600 18.23 0.9515
Ne 614.5 576.6 517.6 613.2 2596 169.6
Ar 3819 3642 3282 3812 1.548 × 104 1093
Be 34.86 30.49 27.17 33.75 158.2 8.952
Mg 1086 1024 920.9 1083 4519 303.0
C 126.0 113.3 101.2 122.8 547.9 33.07
N 203.9 185.6 166.1 200.0 876.7 54.24
O 308.6 284.1 254.6 304.6 1317 83.19
Si 1754 1662 1495 1749 7218 493.9
P 2173 2062 1857 2167 8901 614.5
S 2654 2523 2272 2647 1.083 × 104 753.5
MAPE 7.61% 17.2% 1.14% 331% 72.6%
In other words, Veff (r; [ρ]) + V
P (r; [ρ]) is now the one-
body effective potential for the corresponding Boson sys-
tem with the same electron density. This reduction of
an N -fermion problem to a Boson form is widely imple-
mented in OF-DFT due to its numerical stability and its
easy implementation using existing KS-DFT codes [10].
The associated Pauli potentials for the HQ and LQ
models are immediately obtained by subtraction of the
W KP from Eq. (34) and Eq. (36) respectively. The stan-
dard finite difference method for solving Euler equations
for the TFW models [52] are implemented for the LQ
and HQ models, and the nonlocal terms are evaluated by
Fourier transforms. The choices of radial grids for both of
the real and Fourier space and other detailed numerical
methods are given in Ref. [7]. The kinetic energy for HQ
and LQ models is computed using the Herring pathway
in Eq. (8).
As shown in Table I, the atomic energy calculated by
Figure 3: (Color online) Radial density r2ρ(r) of the Kr atom
using the KS method, the LQ and HQ models, the CAT
model, and the TFλW models (see the inset) with the full
nuclear potential.
 0  0.5  1  1.5  2
 0  0.5  1  1.5  2
TF1/5W
the energy-optimized TF1/5W model is very close to the
KS-DFT, and outperforms all the LR-based models, and
other TFλW models. In Table II, we compare the elec-
tron density at the nucleus ρ(0) for various models. The
TF1/5Wmodel overestimates ρ(0) by about a factor of 4,
while the TFW model underestimates it by about 30%.
The predicted values of ρ(0) for all the LR-based mod-
els are very close to the KS results, and are much better
than the TFλW models.
In Fig. (3), we compare the radial density distribu-
tion r2ρ(r) of the LQ and HQ models to that predicted
by other theories for the Kr atom. Both the TF1/5W
and TFW models predict smooth and structureless ra-
dial density profiles. Using the full Coulombic potential,
all the LR-based models can predict an incipient shell
structure for heavy atoms (Z & 30), and these results
are typical. Since the potential is certainly far beyond
the LR regime, these qualitative results with some sug-
gestion of shell structure are about as good as could be
hoped for. The surprisingly good total energies given in
Table I for the TF1/5W model and the LR-based models
shows that averaged thermodynamic properties are less
sensitive to errors in the KP than is the density profile.
The difference in the results for the LQ, HQ and CAT
models indicates that the LR-based OF theory is being
used outside its range of validity. As shown below, we
gain a significant improvement by using the AILPS to
deal with these difficulties.
B. Ab initio local pseudopotential calculations
As discussed earlier, the use of pseudopotentials in non-
local LR-based OF-DFT can improve the accuracy of the
theory because the weaker pseudopotential is more nearly
in the LR regime, where the theory is designed to be ac-
Figure 4: The smooth target density ρ̃v(r) from Eq. (40), with
parameters given in Table IV for the Si pseudoatom used in
the inverse-KS process, and the valence density ρv(r) pre-
dicted by the LQ and HQ models using the Vps(r) (see Fig.
5) corresponding to ρ̃v(r). The arrow indicates the location
of rc.
 0.01
 0.02
 0.03
 0.04
 0.05
 0.06
 0.07
 0  0.5  1  1.5  2  2.5  3
smooth
Figure 5: The AILPS Vps(r) for Si generated by the target
density ρ̃v(r) in Fig. 4.
 0  0.5  1  1.5  2  2.5  3
curate. Our proposed HQ and LQ models can be used
with any existing AILPS. However, since we want to as-
sess the performance of these models for a wide class of
atomic systems, we describe here a new method for de-
termining reasonable AILPS for general atomic systems.
These pseudopotentials will be used in all our calculations
and can be transfered to other molecular and solid state
environments, but we expect (and find in cases where
comparison can be made) little change if other reason-
able AILPS are used.
Because of the one-to-one mapping between the effec-
tive one-body potential acting on a system of N electrons
and the electron density in the ground-state configura-
tion, it is possible to obtain a unique local one-body po-
tential that generates a given target density ρ̃(r) by using
Figure 6: (Color online) Radial valence density r2ρv(r) of
the Si atom using the KS method and various models using
AILPS. Parameters used for constructing this reference sys-
tem are shown in Table IV. The arrow indicates the location
of rc. Inset: The corresponding radial total density r
2ρ(r),
which is dominated by the core component for r < rc.
 0.05
 0.15
 0.25
 0.5  1  1.5  2  2.5  3
TF1/5W
 0  0.5  1  1.5  2
a KS orbital-based method in an inverse way [23]. To con-
struct an AILPS for a given atom we separate the total
electron density ρ(r) into a “core density” ρ̃c(r), which is
supposed not to vary significantly in other molecular or
solid state environments, and the target“valence density”
ρ̃v(r) where
ρ(r) = ρ̃v(r) + ρ̃c(r). (39)
Because DFT requires only the electron density, we can
take a more general view of what is meant by the core and
valence components than is used in most orbital-based
methods. Here, we directly construct a smooth target
valence density for the Nv = N − Nc valence electrons,
with Nc chosen to be the number of electrons in the noble
gas configuration.
Our proposed target valence density ρ̃v(r) for atoms
equals the full KS density ρKS(r) outside a core of radius
rc, and is designed to be small and slowly varying inside
rc. The functional form we take is
ρ̃v(r) =
tρKS(rc) + a0r
q exp[−rp(a
r2)] r ≤ rc
ρKS(r) r > rc
Figure 4 gives an example of ρ̃v(r) for Si that will be
discussed in more detail below. We find most results are
insensitive to the details of our fitting procedure. Param-
eter values for a variety of atomic systems are given in
the Appendix in Table IV along with the physical and
technical considerations that guided our choice of this
particular form for ρ̃v(r). The Appendix also discusses
some of the general issues that arise in using these atomic
AILPS in other environments.
The local pseudopotential is directly related to the
effective one body potential that reproduces ρ̃v(r) ex-
actly when using the full KS theory. Following previous
work [24, 25], for a given ρ̃v(r), the inverse-KS equations
are solved to get the effective one-body screened poten-
tial Vscr(r). The desired ab initio local pseudopotential
Vps(r) is then obtained by subtracting the Hartree po-
tential and the exchange-correlation potential:
Vps(r) = Vscr(r)− VH(r; [ρ̃v])− Vxc(r; [ρ̃v]). (41)
This relatively expensive procedure to determine Vps(r)
requires the use of orbitals. However it needs to be done
only once for each atom, and the resulting Vps(r) can
then be used in a variety of other environments if the
atomic core densities remain essentially constant.
Once suitable Vps(r) have been determined (by this or
other means) they can be incorporated in Vext(r) in dif-
ferent ways, depending on the particular system of inter-
est. OF-DFT theory can then be used to determine the
valence density ρv(r) in direct analogy to the all-electron
calculations for the full atomic potential in Eq. (38):
▽2 + VH(r; [ρv]) + Vxc(r; [ρv])
+Vext(r) + VP (r; [ρv])}ψv(r) = µψv(r) (42)
where
ρv(r) = |ψv(r)|
2. (43)
A simple and direct test of OF-DFT is to use Eqs. (42)
and (43) for the same atomic system for which Vps(r) was
constructed. Thus we take Vext(r) = Vps(r) for a given
atom as input data in Eq. (42). The valence density
ρv(r) predicted by the OF KPs is determined from Eqs.
(42) and (43), and can be directly compared to the exact
target density ρ̃v(r) for this atomic system given by the
full KS theory.
This is illustrated in Fig. 4, which shows the input
target valence density ρ̃v(r) for Si used in the inverse-
KS process. The total KS density ρKS(r) equals ρ̃v(r)
for r ≥ rc, indicated by the arrow in Fig. 4, and then
increases rapidly for r < rc, reaching a very large value
at the nucleus, ρKS(0) = 1754. In contrast, the proposed
target valence density ρ̃v(r) remains small and relatively
slowly varying inside the core, with ρ̃v(0) = 0.005621.
Also shown in Fig. 4 are the predicted valence densities
ρv(r) for Si given by the LQ and HQ models. Because of
the relatively weak Vps(r) and slowly varying valence den-
sity, both the LQ and HQ KP models predict results very
close to those given by the exact KS treatment of the ki-
netic energy, and perform markedly better than they did
for the all-electron calculations using the full Coulomb
potential. Fig. 5 shows the corresponding AILPS gener-
ated by the inverse KS procedure. It is much smaller in
the core region than the full atomic potential and more
likely to be accurately treated by LR-based methods.
The radial valence density r2ρv(r) of the Si and Ar
atoms predicted by the various methods are shown in
Figs. 6 and 7. The consistency of our OF theory when
pseudopotentials are used is illustrated by the similar-
ity of the density predicted by the HQ and LQ models.
Figure 7: (Color online) Same as in Fig. 6 but for the Ar
atom.
 0.5  1  1.5  2  2.5  3
TF1/5W
 0  0.5  1  1.5  2
The slight deviations from the KS-DFT results for Ar in
the valence region in Fig. 7 are among the largest we en-
countered for all atoms tested, and could be due to the
relatively large number of valence electrons (8) compared
to core electrons (10). Further errors may arise from the
LFWV approximation in Eq. (33).
Once ρv(r) has been determined using OF-DFT, it can
be added to the known input core density ρ̃c(r) to obtain
the predicted total density, since the basic assumption of
our AILPS is that the core density remains unchanged
in different chemical environments. The core density is
defined in Eqs. (39) and (40). The inset in Fig. 6 shows
that the HQ and LQ treatment of the valence density
for the Si atom does not produce noticeable errors in the
total density, as expected, and the shell structure remains
in excellent agreement with the full KS-DFT calculations.
As a further test of OF-DFT we can compare the en-
ergy for the valence density given by the various methods
to the exact valence energy for the target valence density
determined by the inverse KS method. Table III gives
the valence energy values for the KS-DFT, and the LQ
and the HQ models. As can been seen, both the LQ
and HQ models give very good agreement with KS-DFT
and perform significantly better than the other models.
These results show that for this class of relatively weak
pseudopotentials the OF treatment of the KP is quite
satisfactory.
To test of the transferability of the present AILPS, we
also performed calculations for positive ions. Fig. 8 shows
that the ionization energies of various atoms calculated
using KS-DFT and the full atomic densities and those
from the KS-DFT using the valence densities with AILPS
are very similar. Therefore, the present AILPS are quite
transferable to these positive ions.
The ionization energy for models using the AILPS is
obtained by subtracting the valence energies for systems
with Nv and Nv − 1 electrons. Since the core electrons
are assumed to be unaltered in different chemical envi-
Table III: The total valence energy Ev[ρv] using the KS
method, the LQ and HQ models, the CAT model, and the
TFλW models. MAE, the mean absolute error (relative to
the KS method) of various OF models are given at the bot-
tom of their respective columns. Parameters used in Eq. (40)
for such systems are given in Table IV.
KS LQ HQ CAT TF1/5W TFW
Be −0.9914 −0.8955 −0.8950 −0.9583 −1.214 −0.7786
C −6.134 −6.080 −6.100 −6.345 −7.761 −5.266
N −11.04 −11.06 −11.09 −11.32 −13.93 −9.462
O −18.01 −18.09 −18.09 −18.19 −22.44 −15.30
Si −3.771 −3.738 −3.750 −3.869 −4.467 −3.350
P −6.474 −6.432 −6.455 −6.582 −7.385 −5.756
S −10.20 −10.10 −10.14 −10.24 −11.27 −9.023
Ar −21.37 −20.84 −20.91 −20.85 −22.40 −18.56
MAE 0.119 0.103 0.184 1.610 1.312
ronments, the atomic core energy is a constant that can-
cels here or in other similar applications to molecules and
solids. Limitations of the LQ and HQ models are more
evident here, but they do capture the overall periodicity
of the ionization energies well, and perform significantly
better than the other models.
VI. CONCLUSION
In summary, we propose two nonlocal OF KPs that
satisfy exact limits for small and large wavevector per-
turbations and reproduce the exact LR function in the
homogeneous limit. These are the same limits that sev-
eral current KEDFs are designed to satisfy. However,
because of the more local nature of the KP, it is much
easier to satisfy these conditions for the KP than for the
KEDF, and there may be other physical and technical
advantages arising from use the more local KP.
In general, there is no reason to believe that any LR-
based OF-DFT should work well for arbitrary systems
where the model potentials are far beyond the LR regime.
However, most chemical processes involve changes of va-
lence electron densities, which can often be described by
a weak AILPS. Thus the use of a LR-based OF-KP to-
gether with AILPS for atomic systems at least seems well
justified. The small and relatively slowly varying ρv(r)
also provides some justification for our use of the local
FWV kF (r) in Eq. (33).
When the AILPS is used, the valence densities given by
the LQ and HQ KPs are close to those given by the KS
method. Thus the particular integration pathway used
to get the total energy value becomes unimportant. The
simple pathway in Eq. (8) is especially useful, since no
coupling parameter integration is needed.
The proposed models are not only conceptually simple,
but also exact for a model system with a weak poten-
tial and a slowly-varying density. The appearance of the
atomic shell structure was found to be very sensitive to
the accuracy of the proposed KPs. The LR-based LQ and
Figure 8: (Color online) Ionization energies (shown in
Hartree) of the first and the second row atoms using the full
KS method, and various models using AILPS. The mean abso-
lute errors (relative to the full KS method) of various models
using AILPS are: KS (0.1 eV), LQ (1.8 eV), HQ (2.1 eV),
CAT (3.9 eV), TFW (3.1 eV), and TF1/5W (4.3 eV). Ioniza-
tion energies using the TF1/5W model are not shown in the
figure due to its relatively poor performance.
 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18
KSfull
HQ KPs give at best only qualitative indications of shell
structure for full atomic systems, though total energies
are surprisingly good. Still better results for atoms and
ions can be found by focusing on the valence density as
determined by a relatively weak AILPS. While these re-
sults seem promising, improved KPs are needed and fur-
ther investigation is required to see if these ideas can be
usefully applied to other relevant systems like molecules
and solids. Some initial results along these lines will be
reported elsewhere.
Acknowledgments
This work has been supported by the NSF Grant
CHE01-11104, and by the NSF-MRSEC at the University
of Maryland under Grant No. DMR 00-80008. J.D.C.
acknowledges the support from the UMCP Graduate
School Fellowship, the IPST Alexander Fellowship, and
the CHPH Bolck Grant Supplemental Fellowship. We
are grateful to Prof. Emily Carter and members of her
group for many helpful discussions and comments on an
earlier version of this paper.
Table IV: Parameters used in Eq. (A.1) for the target valence
density of various atoms. Here, p = q = 6, and t = 0.1
are used for all systems. The Vps(r) generated from these
parameterized ρ̃v(r) can then be used in OF-DFT.
Nv a0 a1 a2 rc
Li 1 2.983 × 10−4 0.05260 −6.560 × 10−3 2.135
Be 2 0.02655 0.7078 −0.2130 1.370
B 3 0.8032 5.855 −3.531 0.9714
C 4 11.53 30.17 −31.20 0.7429
N 5 98.80 112.8 −180.2 0.5978
O 6 592.7 338.9 −779.8 0.4981
F 7 2750 871.4 −2744 0.4256
Ne 8 1.052 × 104 1993 −8267 0.3707
Na 1 5.234 × 10−5 9.840 × 10−3 −6.610 × 10−4 2.904
Mg 2 9.805 × 10−4 0.04445 −5.002 × 10−3 2.233
Al 3 6.164 × 10−3 0.1273 −0.02056 1.861
Si 4 0.02942 0.3119 −0.06861 1.593
P 5 0.1149 0.6864 −0.1979 1.390
S 6 0.3841 1.390 −0.5100 1.231
Cl 7 1.133 2.632 −1.201 1.103
Ar 8 3.018 4.714 −2.624 0.9985
Appendix: AB INITIO LOCAL
PSEUDOPOTENTIALS
As discussed above, our proposed target valence den-
sity for atoms has the following form:
ρ̃v(r) =
tρKS(rc) + a0r
q exp[−rp(a
r2)] r ≤ rc
ρKS(r) r > rc
(A.1)
Here p and q are taken as even integers. The larger
they are, the smaller and more slowly varying is the va-
lence density near r = 0 but the sharper is the peak
near the core radius rc. As a compromise, we take here
p = q = 6, which generates relatively slowly varying local
pseudopotentials Vps(r). For applications in different en-
vironments, such as molecules or crystals, the core size rc
has to be small to maintain transferability of the atomic
core density. For this reason, we force ρ̃v(r = 0) to be
small by taking a small t. If t = 0, the strict vanish-
ing of the valence density near the nucleus would require
a very repulsive Vps, which is certainly undesirable for
the LR-based OF-DFT. However, if t is too large, there
will exist a long oscillatory tail outside the core in the
corresponding Vps(r). This is an undesirable feature for
transferability to other environments, as will be discussed
below. These two points constrain the value of t and we
use here t = 0.1 for all the atomic systems considered.
The four parameters a0, a1, a2, and rc are determined
by requiring continuity of the function ρ̃v(r) and its first
two derivatives at r = rc, and by satisfying the normal-
ization condition:
Nv = 4π
ρ̃v(r)r
2dr. (A.2)
Here we used the standard noble gas cores to determine
Figure 9: The inverse KS procedure generates very small
oscillations in the tail of rVps(r) for Si (shown here) and
other atoms. The two points where rVps(r) = −4 for Si are
r1 = 2.336 and r2 = 4.576 or r
= 2.139. The arrow indi-
cates the location of r2. To achieve good transferability (see
text), this rVps(r) is modified by setting rVps(r) = −4 when
r ≥ r2 = 4.576. See Fig. 5 for a large scale view.
-4.03
-4.02
-4.01
-3.99
-3.98
-3.97
-3.96
-3.95
 1.6  1.8  2  2.2  2.4  2.6  2.8  3  3.2  3.4
Nv = N −Nc, though other choices could in principle be
made.
To construct our local pseudopotential Vps(r) for
atoms, we first solve the KS equations for an atom with
the full Coulomb potential. With the KS density ρKS(r)
and Nc determined, this construction ensures that as
r → rc, the associated core density smoothly approaches
zero as O(| r − rc |
After generating the parameterized target valence den-
sity ρ̃v(r), the set of inverse-KS equations are solved
to obtain the corresponding one-body screened poten-
tial Vscr(r). In principle, the AILPS Vps(r) is then given
by Eq. (41). Using this, we find with an acceptable Nc
that essentially the same density profiles outside rc are
predicted by the LQ and HQ models for a wide range of
choices of p, q, and t. This would be expected if most fea-
tures of the resulting set of model potentials are within
the different regimes accurately described by the OF KP.
These results clearly show that OF-DFT can give accu-
rate results for this class of physically relevant and rela-
tively weak and slowly varying model potentials.
However, we found that the rVps(r) constructed in this
way for group III to group VIII elements deviates from
−Zv outside the core by a small and long range oscilla-
tion. For the 2s22px atoms, the maximum amplitude of
the oscillation is about 0.18. For the 3s23px atoms, it
ranges from 0.003 for Al to 0.037 for Ar. Similar oscil-
latory tails were observed by Wang and Stott [25], and
are thought to arise from the inability of a local pseu-
dopotential to represent both s and p orbitals of the cor-
responding nonlocal pseudopotentials. These oscillations
are so small in magnitude that they have almost no effect
on the density profile or energy of atoms, but they can
cause transferability problems when used for solids [53].
Following Wang and Stott’s approach, we replace Vps(r)
by −Zv/r at the largest point where rVps(r) = −Zv (see
Fig. (9)). The new Vps(r) can then be used for both
atomic and solid state calculations.
In other chemical environments, such molecules and
solids, the predetermined AILPS Vps(r) centered at each
nucleus are regarded as input data, and the valence den-
sities for other systems can be calculated by OF-DFT.
Results can be checked by using the full KS-DFT. For
example, for molecules and solids, the external potential
Vext(r) is a linear combination of the local atomic pseu-
dopotentials centered at each ion position RI :
Vext(r) =
Vps(r−RI) (A.3)
The valence density can then be determined from Eqs.
(42) and (43), and the total energy of valence electrons
Ev[ρv] = Ts[ρv] + EH [ρv] + Exc[ρv] +
ρv(r)Vext(r)dr
(A.4)
For the proposed LQ and HQ KPs, the value of Ts[ρv]
in Eq. (A.4) is determined by Eqs. (8) or (11) with ρ(r)
replaced by ρv(r).
[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864
(1964).
[2] R. G. Parr and W. Yang, Density-Functional Theory of
Atoms and Molecules, (Oxford University Press, 1989).
[3] R. M. Dreizler and E. K. U. Gross, Density Functional
Theory: An Approach to the Quantum Many Body Prob-
lem, (Springer-Verlag, Berlin, 1990).
[4] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[5] L. J. Sham and W. Kohn, Phys. Rev. 145, 561 (1966).
[6] R. A. King and N. C. Handy, Phys. Chem. Chem. Phys.
2, 5049 (2000); R. A. King and N. C. Handy, Mol. Phys.
99, 1005 (2001).
[7] J.-D. Chai and J.D. Weeks, J. Phys. Chem B 108, 6870
(2004).
[8] Y. A. Wang, N. Govind, and E. A. Carter, Phys. Rev. B
58, 13465 (1998).
[9] Y. A. Wang, N. Govind, and E. A. Carter, Phys. Rev. B
60, 16350 (1999).
[10] See e.g., Y. A. Wang and E. A. Carter, in Theoreti-
cal Methods in Condensed Phase Chemistry , edited by
S. D. Schwartz, Progress in Theoretical Chemistry and
Physics, (Kluwer, Boston, 2000) p. 117, and references
therein.
[11] B. Zhou, V. L. Ligneres, and E. A. Carter, J. Chem.
Phys. 122, 044103 (2005).
[12] E. Chacón, J. E. Alvarellos, and P. Tarazona, Phys. Rev.
B 32, 7868 (1985).
[13] P. Garćıa-González, J. E. Alvarellos, and E. Chacón,
Phys. Rev. A 54, 1897 (1996).
[14] P. Garćıa-González, J. E. Alvarellos, and E. Chacón,
Phys. Rev. A 57, 4857 (1998).
[15] M. Pearson, E. Smargiassi, and P. A. Madden, J. Physics:
Condens. Matter 5, 3221 (1993); E. Smargiassi and P. A.
Madden, Phys. Rev. B 49, 5220 (1994); M. Foley and P.
A. Madden, Phys. Rev. B 53, 10589 (1996).
[16] F. Perrot, J. Phys.: Condens. Matter 6, 431 (1994).
[17] J. A. Alonso and L. A. Girifalco, Phys. Rev. B 17, 3735
(1978); M. D. Glossman, L. C. Balbás, and J. A. Alonso,
Chem. Phys. 196, 455 (1995).
[18] C. Herring, Phys. Rev. A 34, 2614 (1986).
[19] L.-W. Wang and M. P. Teter, Phys. Rev. B 45, 13196
(1992).
[20] N. Choly and E. Kaxiras, Solid State Commun. 121, 281
(2002).
[21] L. H. Thomas, Proc. Cambridge Phil. Soc. 23, 542
(1927); E. Fermi, Z. Phys. 48, 73 (1928).
[22] See e.g., R. Evans, in Fundamentals of Inhomogeneous
Fluids, edited by D. Henderson (Dekker, New York,
1992), p. 85.
[23] Y. Wang and R. G. Parr, Phys. Rev.A 47, R1591 (1993).
[24] B. Zhou, Y.A. Wang, and E.A. Carter, Phys. Rev. B 69,
125109 (2004).
[25] B. Wang and M. J. Stott, Phys. Rev. B 68, 195102
(2003).
[26] A. Holas and N. H. March, Phys. Rev. A 66, 066501
(2002); I. Lindgren and S. Salomonson, Phys. Rev. A
67, 56501 (2003).
[27] L. R. Pratt, G. G. Hoffman, and R. A. Harris, J. Chem.
Phys. 88, 1818 (1988); L. R. Pratt, G. G. Hoffman, and
R. A. Harris, J. Chem. Phys. 92, 6687 (1990); G. G.
Hoffman and L. R. Pratt, Mol. Phys. 82, 245 (1994).
[28] Y.-G. Chen and J. D. Weeks, J. Chem. Phys. 118, 7944
(2003)
[29] J. P. Perdew and S. Kurth, in Density Functionals: The-
ory and Applications, edited by D. Joubert, Lecture
Notes in Physics, (Springer-Verlag, Berlin, 1998) p. 8,
and references therein.
[30] L. J. Sham, Phys Rev. A 1, 969 (1970).
[31] M. Levy and J. P. Perdew, Phys. Rev.A 32, 2010 (1985).
[32] F. W. Averill and G. S. Painter, Phys. Rev. B 24, 6795
(1981).
[33] W. Yang, Phys Rev. A 34, 4575 (1986).
[34] C. F. von Weizsäcker, Z. Physik 96, 431 (1935).
[35] J. Lindhard, K. Dan. Vidensk. Selsk. Mat. Fys. Medd.
28, 8 (1954).
[36] W. Jones and W. H. Young, J. Phys. C 4, 1322 (1971).
[37] Y. Tomishima and K. Yonei, J. Phys. Soc. Jpn. 21, 142
(1966).
[38] E. H. Lieb, Rev. Mod. Phys. 53, 603 (1981).
[39] G. K.-L. Chan, A. J. Cohen, and N. C. Handy, J. Chem.
Phys. 114, 631 (2001).
[40] P. K. Acharya, L. J. Bartolotti, S. B. Sears, and R. G.
Parr, Proc. Natl. Acad. Sci. USA 77, 6978 (1980).
[41] J. L. Gázquez and J. Robles, J. Chem. Phys. 76, 1467
(1982).
[42] L. J. Bartolotti and P. K. Acharya, J. Chem. Phys. 77,
4576 (1982).
[43] P. K. Acharya, J. Chem. Phys. 78, 2101 (1983).
[44] D. J. Singh, Phys. Rev. B 48, 14099 (1993).
[45] J.-D. Chai, Ph.D. dissertation, University of Maryland,
2005.
[46] Instead of solving the second-order differential equation
numerically for the corresponding weight function ω(x)
(see Eq. (14) of Ref. [13]), we used the parametrized form
in the appendix of Ref. [13]. Although the parametriza-
tion was made for β = −1/2 (note the sign change of β
in Ref. [13]), the β dependence of ω(x) is small for small
β, as argued in Ref. [14]. Our results for the CAT model
using this parametrized ω(x) are very close to their re-
ported results.
[47] P. A. M. Dirac, Proc. Cambridge Phil. Soc. 26, 376
(1930).
[48] D. M. Ceperley, Phys. Rev. B 18, 3126 (1978); D. M.
Ceperley and B.J. Alder, Phys. Rev. Lett. 45 , 566
(1980).
[49] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048
(1981).
[50] N. H. March, Phys. Lett. A 113, 476 (1986).
[51] A. Holas and N. H. March, Phys. Rev.A 44, 5521 (1991).
[52] A. M. Abrahams and S. L. Shapiro, Phys. Rev. A 42,
2530 (1990).
[53] D. M. Bylander and L. Kleinman, Phys. Rev. B 55, 9432
(1997).
