Eur. Phys. J. B (2018) 91: 189
https://doi.org/10.1140/epjb/e2018-90111-0

THE EUROPEAN
PHYSICAL JOURNAL B

Regular Article

On-top density in the nonlinear metallic screening and its
implication on the exchange-correlation energy functional?
Yasutami Takada aInstitute for Solid State Physics, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8581, Japan

Received 1 March 2018 / Received in final form 8 June 2018
Published online 29 August 2018c EDP Sciences / Società Italiana di Fisica / Springer-Verlag GmbH Germany, part of Springer Nature,
2018
Abstract. In comparison with the accurate data on the on-top electron density n(0) in the protonembedded electron gas with the density parameter rs in the range 1−12 obtained by diffusion Monte
Carlo (DMC) simulations, we have successfully constructed an alternative form of the exchange-correlationenergy functional in the density functional theory by imposing the constraint due to the cusp theoremon the well-known Perdew–Burke–Ernzerhof (PBE) functional. Although PBE does not, our functional,
referred to as the cusp-corrected PBE (ccPBE), reproduces the DMC data on n(0) in the entire range of rs .

1 Introduction
An atom, especially hydrogen, immersed into the otherwise homogeneous electron gas (EG) has been investigatedfor more than four decades not only in the density functional theory (DFT), mostly in its local-density approximation (LDA) [1–14], but also in various forms of manybody theories [15–21], including diffusion Monte Carlo(DMC) [22] and variational Monte Carlo (VMC) [23] simulations. The primary motivation of those studies is toconstruct an appropriate theory for the nonlinear responseof metallic electrons to an impurity point charge +Ze,
but the basic physical concept with which they were concerned remains the same as that in the linear-response
theory, known as Thomas–Fermi (TF) [24,25] (or Debye
and Hückel [26]) screening of the impurity charge with ashort screening length λTF which is about the same as
kF−1 , where kF is the Fermi wave number of EG.
Recently, by studying a proton (the case of Z = 1)
embedded in EG with use of both LDA and DMC, thepresent author has gained a new insight into this problem [27]; the concept of Kondo screening of the spin ofhydrogen with a long screening length λK (kF−1 ) [28] isfound to be relevant to this system and a sharp transition from TF to Kondo screening is shown to exist withthe decrease of the metallic electron density n0 (=kF3 /3π 2 )from the high-density limit. At the same time, the resultsin DMC are found to be well approximated by thosein LDA in the density region characterized by Kondo
screening because of the slowly-varying nature of the?
Contribution to the Topical Issue “Special issue in honorof Hardy Gross”, edited by C.A. Ullrich, F.M.S. Nogueira,A. Rubio, and M.A.L. Marques.
ae-mail: takada@issp.u-tokyo.ac.jp

electron density distribution n(r) around the proton due
to the long λK .
In the high-density region characterized by TF screening, on the other hand, a relatively large difference canbe seen in n(r) between DMC and LDA. In particular,
an unexpected feature of n(r) is found at the proton
position or the on-top density n(0); in DMC, nDMC (0),is lower than that in LDA, nLDA (0), at high densities,
namely, rs < 1.66 with the conventional density parameter rs ≡ (αkF aB )−1 , while the opposite is the case forrs > 1.66. Here we define α = (4/9π)1/3 ≈ 0.5211 and aB
is the Bohr radius. (We will use atomic units hereafter.)
According to a physical argument [29], we obtain largern(0) for stronger exchange-correlation (xc) effect, implying that as long as we believe that nDMC (0) is sufficiently
accurate, LDA is found to provide a too strong xc effect for
rs < 1.66 but a too weak one for rs > 1.66. This interesting crossover behavior with the increase of rs in describing
the xc effect in LDA has never been known.
Because n(r) varies very weakly even for r < λTF in
densely packed systems such as the high-density EG for
rs < 1.66, it is natural to expect that a small densitygradient correction to LDA will be enough to obtaina result of n(0) much better than LDA, but actually
the situation becomes worse in the generalized gradientapproximation (GGA) in the Perdew–Burke–Ernzerhof(PBE) version [30,31]; namely, for rs < 1.66, the difference of PBE from DMC becomes larger than that of LDA.
One might imagine that not PBE but the accurate gradient expansion [32] as included in PBEsol [33,34] is needed
to obtain better n(0), but this is not the case; no improvement on LDA is achieved even in PBEsol. Thus we cometo notice that it is a nontrivial work to reproduce nDMC (0)for the case of rs < 1.66 in the framework of GGA.

Page 2 of 9

Eur. Phys. J. B (2018) 91: 189

In pursuit of a key ingredient to improve on PBE in
the present problem with retaining exact conditions whichmake PBE reliable, as listed in Table 1 in reference [35], wecome across the importance of the cusp theorem [36–38]
which dictates that n(r) near the impurity atom behaves
rigorously in such a manner asn(r) −−−→ ncusp (r) ≡ n(0) exp(−2Zr),
r≈0

(1)

so as to make a compromise with the singular Coulomb
potential term −Z/|r|. Although nLDA (r) satisfies
equation (1), n(r) in PBE or PBEsol does not, indicatingthat the worse performance of PBE/PBEsol in determining n(0) might originate from the violation of the cusp
theorem.
Generally it is not believed that we can make the cusp
theorem obeyed in the framework of GGA [35,39] and it
is usually thought that some form of meta-GGA [40–42]
is needed to satisfy it. Therefore the inclusion of the cusp
theorem into a GGA-based scheme is really a challenge. In
this paper we take up this challenge and set the goal of thispaper in the following way; we just try to modify the spinresolved xc energy functional E xc [nσ ] in PBE by imposing
the constraint due to the cusp theorem in addition to theexact conditions already obeyed by PBE and then we tune
up some free parameters involved in the modified E xc [nσ ]so as to reproduce nDMC (0) in the wide range of rs , i.e.,
1 ≤ rs ≤ 12 where the DMC data are available.
We will leave a comprehensive test of this modified
E xc [nσ ] (which will be referred to as “cusp-corrected”
PBE or ccPBE) for a variety of real materials for the
future, but because ccPBE provides the different resultsof n(0) from those in PBE only for rs < 1.66, ccPBE andPBE will give, more or less, similar results for almost
all real materials. One important exception is the solidhydrogen under very high pressures [43,44] in which
1.1 < rs < 1.7. Thus ccPBE may be expected to be useful
only for solid hydrogen.
This paper is organized as follows: In Section 2, we
introduce the system to be treated, explain the calculation methods, and account for the issues arisen from thedata calculated on n(0). In Section 3, we construct ccPBE
and give the calculated results for n(0) in ccPBE in comparison with those in DMC. Finally in Section 4, we givea summary of this paper and make several comments.

The Hamiltonian H for electrons in this system is given as
H =−

X ∇2
i

i

2

+

X
1X
1
+
vext (ri ),2|ri −rj |
i

(2)

i6=j

where the external potential working on an electron
vext (r) is composed of the potential from the nucleus andthat from the positive background, written asZ
N − Z 3R2 − r2
θ(R − r)
−
r
2
R3
N −Z
−
θ(r − R),rvext (r) = −(3)

with r = |r| and θ(x) the Heaviside step function. In
solving equation (2), we impose the fixed boundary
condition to make the wave function vanish at |ri | = R.
From a computational point of view, this boundary
condition is indispensable to obtain rapidly and stably
convergent results in the closed-shell condition.
2.2 DFT and the Kohn–Sham scheme
In DFT, the spin-resolved ground-state density nσ (r) for
H in equation (2) is rigorously determined by the map to
a noninteracting reference system which is solved by the
Kohn–Sham (KS) equation, written as−∇2 /2 + vσKS (r) φiσ (r) = εiσ φiσ (r),(4)

where εiσ and φiσ are the energy level and the normalized
wave function for KS orbital i and spin σ, respectively,
and vσKS (r) is the KS potential, given byvσKS (r) = vext (r)+

Z

dr 0

n(r 0 )
+v xc (r; [nσ ]),|r−r 0 | σ(5)where vσxc (r; [nσ ]) is derived from E xc [nσ ] through thefunctional derivative as
vσxc (r; [nσ ]) = δE xc [nσ ]/δnσ (r).

(6)

With use of the lowest-Nσ KS orbitals, nσ (r) is given bynσ (r) =Nσ
X|φiσ (r)|2 ,

(7)

i=1

2 Atom embedded in the jellium sphere
2.1 Hamiltonian
Because Monte Carlo simulations can treat only a finite
number of electrons, let us consider not bulk jellium but
a jellium sphere of radius R and average density n0 andthen put a neutral atom of atomic number Z at r = 0(the center of the sphere). The number of electrons in
the jellium sphere is 4πR3 n0 /3 = (R/rs )3 , so that the
total electron number N is equal to Z+(R/rs )3 , satisfying
global neutrality, from which we obtain R = (N −Z)1/3 rs .

and n(r) is the sum of n↑ (r) and n↓ (r). The
P spindensity nσ (r) and consequently Nσ with N = σ Nσshould be determined by the self-consistent solution ofequations (4)–(7), together with the fixed boundary condition
φiσ (r) = 0,

(8)

at r = R = (N − Z)1/3 rs . This boundary condition is
imposed to make a direct comparison of the results in
DFT-based schemes with those in DMC.

Eur. Phys. J. B (2018) 91: 189

Page 3 of 9

where u = u(r) and v = v(r) are defined, respectively, as

2.3 LSDA
In order to implement the above KS scheme, we need toknow some concrete form of E xc [nσ (r)]. In the local-spin
density approximation (LSDA), it is given byZxc

E [nσ (r)] =
rs (r), ζ(r) ,dr n(r)unifxc(9)

where n(r) = n↑ (r) + n↓ (r) and unifxc (rs , ζ) is the xc
energy per electron for the homogeneous electron gas withthe density parameter rs = (3/4πn)1/3 and the spin polarization ζ = (n↑ −n↓ )/n. Usually, unifxc (rs , ζ) is divided into
two parts; the exchange part unifx (rs , ζ) and the correlation part unif(rs , ζ), both of which are concretely givencin reference [45], but we can simply write unifx (rs , ζ) as
4/3

unif
unif
x (rs , ζ) = x (rs )(1 + ζ)+ (1 − ζ)
2

4/3

,

(10)

2/3
/rs .
with unifx (rs ) = −(3/4)(3/2π)

In GGA, E xc [nσ (r)] is given as a functional of not onlynσ (r) but also its first derivative ∇nσ (r). In its PBE
version, E xc [nσ (r)] is assumed to be
E x [2n↑ (r)]+E x [2n↓ (r)]+E c [nσ (r)], (11)2with the exchange energy functional E x [n(r)], written as
E x [n(r)] =Z

dr n(r)unifx (rs )Fx (s),(12)

where s = s(r) is the normalized derivative, defined bys(r) =|∇n(r)||∇n(r)|=,
2kF (r)n(r)2(3π 2 )1/3 n(r)4/3(13)

with kF (r) = [3π 2 n(r)]1/3 and Fx (s) is given by
Fx (s) = 1 + κ −κ
1 + µPBE s2 /κ∇n(r)·∇|∇n(r)|∇2 n(r), v(r) =. (16)
3
2[2kF (r)] n(r)[2kF (r)]2 n(r)On the other hand, the correlation energy functional
E c [nσ (r)] in equation (11) is given byZEc [n(r)] = dr unif(rs , ζ) + H(rs , ζ, t) ,(17)cwhere the functional H(rs , ζ, t) is defined as



βMB 2
1+At2
H(rs , ζ, t) = γφ ln 1+, (18)tγ
1+At2 +A2 t4
3√
with introducing t = (3π 2 /16)1/3 s/ rs φ(ζ) and the function φ(ζ) defined as φ(ζ) = [(1 + ζ)1/3 + (1 − ζ)1/3 ]/2. HereβMB = 0.066725 [46] and the functional A is given asA=

2.4 PBEE xc [nσ (r)] =u(r) =,(14)

with κ = 0.804 and µPBE = 0.21951. By using E x [n(r)] in
equation (12), we can derive vσx (r) the exchange part ofvσxc (r; [nσ ]) for a spin-σ electron as
δE x [n(r)]δn(r) n(r)=2nσ (r)4∂Fx (s)= unif(r)
Fx (s) − v
s
x
3
s∂s




4
∂ ∂Fx (s)− u − s3
,
3
∂s
s∂s

vσx (r) =

(15)

βMB 
−1
exp[−unif(rs , ζ)/γφ3 ] − 1
,
c
γ

(19)

with γ = (1 − ln 2)/π 2 . By the functional derivative of
E c [nσ (r)] with respect to nσ (r), we obtain vσc (r). The
concrete form for vσc (r) is suppressed here.
In PBEsol, the same forms for the exchange and correlation energy functionals are adopted with the replacement of µPBE and βMB by µGE (=10/81) [32] and 0.046,
respectively.
2.5 DMCThe detailed account of the procedure for DMC is given
in reference [27] and thus we will not recapitulate it here,
but the point is that the only approximation involvedin DMC is the so-called “fixed-node approximation”. As
explained in reference [27], we consider that unphysical node-position dependent effects will be removed byextracting the N -independent results, because the node
positions depend on N in the fixed boundary condition.
As for the on-top density n(0), the N -independent
results are found to be obtained, if N becomes as large as
about 60 for the proton-embedded EG. Furthermore, thecenter of the sphere is very much separated from any postulated node positions, implying that n(0) is the physicalquantity least affected by the fixed-node approximation.
For those reasons, it is well expected that DMC provides
accurate, if not exact, results for n(0).
2.6 On-top density in LDA, PBE, and DMC
We have applied DMC to the system described by H inequation (2) with Z = 1 and, as reported in reference [27],
we have obtained convergent results at N = 58 for rs ≤ 2.6(the TF-screening region) and N = 60 otherwise (theKondo-screening region). In the stably convergent closedshell condition, the doubly-degenerate 3s energy level
corresponding to the Kondo singlet state is situated justabove (below) the Fermi level for low- (high-)rs systems,

Page 4 of 9

Eur. Phys. J. B (2018) 91: 189

Table 1. On-top density n(0) in atomic units for the
proton-embedded electron-gas sphere with the total electron number N = 58 for rs = 1.0−2.6 and N = 60
otherwise.

Fig. 1. (a) Density distribution n(r) calculated in LDA,
PBE, and DMC for the proton-embedded electron-gas spherewith the total electron number N = 58 for rs = 1 and 2,
and N = 60 otherwise. (b) The corresponding normalized
derivative function s(r) in PBE.

leading to the difference in N by 2 between the case ofrs ≤ 2.6 and that of rs > 2.6. In order to make a directcomparison with those DMC results, both LSDA and PBE
have been performed in exactly the same situation as for
N and the boundary condition at each rs . Note that there
is no difference between LSDA and LDA in the present
system, because the ground states are always found to beparamagnetic. Thus we will simply write “LDA” hereafter,even though the actual calculations are done in LSDA.
The obtained results of n(r) in both LDA and PBE
are in good agreement with those in DMC, especially forr larger than 2aB , as seen, for example, in Figure 5 inreference [27], irrespective of either TF or Kondo region.
Relatively speaking, for r less than about aB , however,there are rather large differences among DMC, LDA, andPBE, as shown explicitly in Figure 1a and the largest
deviation occurs at r = 0. Thus it is important to makea detailed quantitative comparison of the on-top density
n(0) among those calculation methods in order to assess
the performance of DFT-based schemes in reference to
DMC.
In view of equation (1), n(r) changes linearly with
the increase of r in semi-log plot, as long as r is lessthan about 0.3aB . (In dense systems like rs less thanabout 2, this critical value for the cusp theorem rcusp
becomes smaller; it may be safe to take 0.1aB for rcusp
at rs = 1.) This linear property in semi-log plots isvery useful in estimating n(0) in DMC. The results so
obtained for n(0) in each scheme are given in Table 1, fromwhich we find that for rs less than about 1.66 (RegionI), nDMC (0) < nLDA (0) < nPBE (0), while for larger rs(Region II which includes the TF-Kondo transition point),
nDMC (0) > nPBE (0) > nLDA (0). Note that this interesting

rs

LDA

PBE

DMC

1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.6
2.7
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0

0.91886
0.69298
0.56659
0.489900.44076
0.40806
0.385700.359020.45749
0.425020.369220.34656
0.33550
0.329390.32573
0.323390.321820.32073
0.31995

0.92150
0.697020.572170.49706
0.44946
0.41815
0.39698
0.37198
0.468230.43576
0.37950
0.356390.34551
0.34027
0.337960.33718
0.337110.337310.337590.894 ± 0.034
0.674 ± 0.023
0.556 ± 0.017
0.486 ± 0.019
0.450 ± 0.014
0.418 ± 0.010
0.402 ± 0.009
0.382 ± 0.008
0.476 ± 0.007
0.442 ± 0.006
0.388 ± 0.006
0.359 ± 0.005
0.349 ± 0.003
0.347 ± 0.003
0.345 ± 0.002
0.343 ± 0.002
0.342 ± 0.001
0.341 ± 0.001
0.339 ± 0.001

crossover point from Region I to Region II is situatedin the density region in which the solid hydrogen andrelated materials under high pressures are involved, i.e.,
1.1 < rs < 1.7 [43,44], making the present assessment
relevant and important in studying physics of the solidhydrogen in the framework of DFT.
Incidentally, in Figure 1b, the results for the normalized derivative s(r) in PBE defined in equation (13),corresponding to those of n(r) in Figure 1a, are plotted,
revealing the interesting fact that in Region I, s(r) always
stays less than 0.43, but in Region II, it becomes larger
than that value. It must also be noted that in the Kondoscreening regime in Region II, the behavior of s(r) is muchdifferent from that in the TF-screening regime, providing another piece of evidence for the qualitative difference
between those two regimes of screening.

3 Proposal of cusp-corrected PBE
3.1 Violation of the cusp theorem in PBE
Confronted with the interesting behavior of the difference between PBE and DMC with the increase of rs in
Table 1, we have made various trials to construct a new
xc energy functional in GGA so that n(0) in DMC canbe well reproduced in the entire range of rs , mostly byjust modifying Fx (s) from the original one in PBE, as is
usually the case in most other modifications of E xc [nσ (r)]from PBE, such as WC [47–49]. Incidentally, there is no
problem in Region II; LDA already provides reasonably
good n(0) and PBE improves much on it, but it is by nomeans easy to obtain n(0) in similar accuracy in Region I.
Thus we will focus on that region in the following.
In Region I, s(r) is less than 0.43 and thus we needsome new insight into E xc [n(r)] in this small-s range. In

Eur. Phys. J. B (2018) 91: 189

Page 5 of 9

pursuit of the new ingredient needed for improving on
the PBE energy functional, we have paid attention to thecusp theorem; as mentioned in Section 1, the cusp behavior in equation (1) is correctly reproduced in LDA, but it
is usually not the case in GGA due to the appearance of asingular term −δZ/|r| in the exchange-correlation potential v xc (r) near the nucleus, in addition to the external
singular term −Z/|r|. In the presence of this additional
singular term, the cusp behavior is not determined by Z
but Z + δZ, leading to the relative error in proportion toδZ/Z.
With the use of equations (1), (13), and (16), we findthat for r ≈ 0, s(r), u(r), and v(r) behave, respectively,ass(r) ≈ sc ,u(r) ≈ s3c ,v(r) ≈ s2c −

s2c
,
Zr

(20)

with sc ≡ Z/[3π 2 n(0)]1/3 . Thus for r ≈ 0, the singularcontribution to vσx (r) in equation (15) comes only from theterm in proportion to v(r). More explicitly, the singularterm can be written as
3 ∂Fx (sc ) 1
s2c ∂Fx (sc )unif c, (21)=−−x (rs ) −Zr
sc ∂sc4π ∂scrwith rsc ≡ [3/4πn(0)]1/3 . Similarly, the singular term invσc (r) is written in the form of equation (21) with thereplacement of Fx (sc ) by Fc (sc ), defined as
Fc (sc ) =unif(rsc , 0) + H(rsc , 0, tc )
c
c
unif
x (rs )

(22)√
with tc = (3π 2 /16)1/3 sc / rsc . Then δZ is given by
δZ =

3
3 ∂
(Fx +Fc ) =
s
4π ∂s
2π



∂Fx
∂Fc
+ 2
2
∂s
∂s


,

(23)

evaluated at the cusp position r = 0 with rs = rsc , s = sc ,and ζ = 0.
For the case of Z  1, n(0) is well approximated byeither Z 3 /π in the strong-correlation limit or 2Z 3 /π in the
weak-correlation limit. Then we obtain sc and rsc , respectively, as either (3π)−1/3 ≈ 0.473 and (3/4)1/3 Z −1 or(6π)−1/3 ≈ 0.376 and (3/8)1/3 Z −1 in each limit, implying that sc is in the range (0.376, 0.473) and rsc  1.
However, not only in the present atom-embedded EG but
also in atoms, molecules, and solids in which the conditionof Z  1 is not always satisfied, sc varies in the range from
0.32 to 0.473, still a relatively small range of s around 0.4.
Now, let us take Fc as FcPBE the one given in PBE. Thenthe second component in equation (23) or the function∂FcPBE /∂s2c is concretely known as a function of sc withrsc set equal to (9π/4)1/3 sc /Z for each Z. In Figure 2, thisfunction (or actually its negative, −∂FcPBE /∂s2c ) is plotted
as a function of sc , from which we see that if sc were
zero (or at least very small), the cusp condition would be(almost) fulfilled in PBE, because ∂FcPBE /∂s2c (which isequal to −µPBE at sc = 0, irrespective of Z) is cancelledby ∂FxPBE /∂s2c = µPBE /(1+µPBE s2c /κ)2 . In fact, in theFig. 2. “Target function” or the partial derivative of −Fc
with respect to s2c in PBE plotted as a function of sc for various Z with rs = (9π/4)1/3 sc /Z. For comparison, ∂Fx /∂s2c in
PBE is also plotted by the dotted-dashed curve. In actual cusp
positions, the value of sc is in the range 0.32−0.473.

original PBE, µPBE is so determined as to satisfy thiscondition at sc = 0, a relation to intimately connect FxPBE
with FcPBE . In the actual cusp region in which sc is about
0.4, however, δZ is not small enough and thus the cusp
theorem is violated in PBE; the relative error δZ/Z isabout 1.8% and 1.2% for Z = 1 and 2, respectively, andless than 1% for Z ≥ 3.
The above observation inspires us that if we can modify
Fx (s) so as to cancel ∂Fx (s)/∂s2 with ∂FcPBE /∂s2 for s
not at a single point of s = 0 but in the entire cusp regionof (0.32, 0.473), we can always make δZ vanish at the cusp
point (and thus impose the cusp theorem), whatever value
for sc is determined in the self-consistent calculation ofn(r). This constitutes the main idea of this paper. Of
course, because −∂FcPBE /∂s2 , which will be called “thetarget function” hereafter, depends on Z, we have to treat
its Z dependence appropriately, but for the time being, wetake the target function at Z = 1. Then, for other valuesof Z the cusp theorem will be violated, but in this casethe relative error δZ/Z becomes much smaller than thatin PBE; at the most, it is about 0.36% for Z = 2 or 3.
A formally better scheme to impose the constraint due
to the cusp theorem for any Z will be mentioned in
Section 4. As for the choice of Fc , we have examined thecase of Fc with βMB in equations (18) and (19) replaced byeither 0.046 as in PBEsol or the more refined rs -dependent
one, β(rs ), expressed as [41,50]β(rs ) = βMB

1 + 0.1rs,
1 + 0.1778rs(24)

but we find that no appreciable difference is seen in the
final results for n(0).
3.2 Exchange energy functional in ccPBE
In order to construct Fx (s) in accordance with the abovementioned idea to fulfill the cusp theorem under theassumption that the correlation energy functional is set
equal to FcPBE , we have examined a variety of possible

Page 6 of 9

Eur. Phys. J. B (2018) 91: 189

forms to arrive at the following Fx (s) which is given as
the sum of three terms:
Fx (s) = F0 (s) + F1 (s) + F2 (s),(25)

where F0 (s) is basically the one only slightly modified
from the original form in PBE as
F0 (s) = A0 + A1 −A1,
1 + µ(p) p/A1

(26)

where p ≡ s2 and µ(p) is assumed to beµ(p) = µ1 + (µ0 − µ1 ) exp(−p/s20 ).

(27)

The function F1 (s) is assumed to be
F1 (s) = B0 exp(−p2 /s41 ),(28)

in order to satisfy the exact gradient expansion (GE) of
Fx (s) in the limit of s → 0, known as [51]Fx = 1+µGE p+

146 2 73v −p v+Dp2 +O(∇6 ),
2015
405

(29)

where v is defined in equation (16) and the coefficient D
vanishes according to the best numerical estimate. The
function F2 (s) is so introduced as to impose the constraint due to the cusp theorem; namely, ∂Fx (s)/∂s2 is
set equal to −∂FcPBE /∂s2 with Z = 1 for s in the range(0.32, 0.473). The actual procedure is to begin with theassumption of F2 (s) in the form of
"

6
X
p2
F2 (s) = 4 C0 +
Ci
s2
i=1



p
s22

i #

A0
A1
µ0
µ1
s0
B0
s1
C0
C1
C2
C3
C4
C5
C6
s2

Z=1

Z=2

Z=3

1.036
0.768
0.12345679
0.13170898
1.20–0.0360.180–0.006933655–0.011363996
0.010829969–0.003780562
0.000643956–0.0000554048
0.00000197693
0.142

1.036
0.768
0.12345679
0.13170898
1.20–0.0360.180–0.007332614–0.036651747
0.039997327–0.014587676
0.002663050–0.0002416611
0.00000954069
0.144

1.036
0.768
0.12345679
0.13170898
1.20–0.0360.180–0.007538429–0.038238364
0.046200149–0.017395384
0.003272557–0.0003038872
0.00001237041
0.145

Then, v is approximately given by v = v0 s2 with a coefficient v0 which is calculated by taking the average of
1/Zr by the weight of n(0) exp(−2Zr) in the range of0 ≤ r ≤ rcusp with rcusp which is the critical value for r
satisfying the cusp condition, as introduced in Section 2.6.
More specifically, v0 is calculated asR 
dr exp(−2Zr)/Zr
1=1− R
v0 = 1 −Zr
dr exp(−2Zr)=−exp(−p/s22 ).

(30)Then, under given values for C0 and s2 , we determine
the six coefficients, C1 , . . . , C6 , so as to satisfy the abovementioned condition for fulfilling the cusp theorem.
There are still nine parameters, s0 , s1 , s2 , A0 , A1 , µ0 ,µ1 , B0 , and C0 , to be fixed, but they cannot be chosen independently; there are four important constraints;
in the limit of s → ∞, there is the Lieb–Oxford upperbound 1 + κ with κ = 0.804 for Fx (s) [52], leading to thecondition oflim Fx (s) = A0 + A1 = 1 + κ.

s→∞

Table 2. Set of parameters to specify Fx (s) in ccPBE. In
order to satisfy the cusp condition, three target functions
corresponding to Z = 1, 2, and 3 are considered. Note that
only the parameters C0 , C1 , . . . , C6 , and s2 depend on Z.

(31)

In the limit of s → 0, we should respect equation (29),but because the functional Fx in GGA is assumed to bea function of a single variable s, we need to derive anapproximate expression for v in terms of s in order tomake use of equation (29). As in equation (20), by theuse of the definitions of s and v in equations (13) and(16), respectively, and the behavior of n(r) in equation (1)near the nucleus at which the electron density varies most
rapidly, we obtain


Z2
1
Z
and v = 2 1 −
.
(32)
s=
kF
kF
Zr

ρ20 /2,exp(ρ0 ) − 1 − ρ0 − ρ20 /2(33)

with ρ0 = 2Zrcusp . Because rcusp is about 0.1 or larger andthe case of Z = 1 is considered here, we take ρ0 as 0.20
tentatively in the following. Then, we obtain the small-s
expansion of Fx (s) in the following way:Fx (s) = 1 + µGE p + µ4 p2 ,(34)

with µ4 = (146/2015)v02 − (73/405)v0 ≈ 17.22612.
In accordance with the small-p expansion inequation (34), terms at each order, O(p0 ), O(p), orO(p2 ), should satisfy
Fx (0) = 1 = A0 + B0 ,∂Fx (0)= µGE = µ0 ,∂p
1 ∂ 2 Fx (0)µ1 − µ0
µ20
B0
C0
=
µ
=
−
− 4 + 4,
4
2
2 ∂2p
s0
A1
s1
s2

(35)
(36)
(37)

respectively. By use of equations (31), (35)–(37), theparameters, A0 , A1 , µ0 , and C0 can be determined under
given values for the rest of the parameters.
By comparing the calculated results for the singleproton embedded electron-gas sphere in ccPBE withthose in DMC, we can determine an appropriate setof parameters providing sufficiently good results. The

Eur. Phys. J. B (2018) 91: 189

Page 7 of 9

Table 3. On-top density n(0) in ccPBE for the protonembedded electron-gas sphere with the total electron
number N = 58 for rs = 1.0−2.6 and N = 60 otherwise.
The parameter sets are used for the target function witheither Z = 1 or Z = 2. The errors are given as the relative
ones in % with respect to the DMC results.
rsFig. 3. Fx (s) in ccPBE determined in reference to three different target functions with Z = 1, 2, and 3. For comparison,
we also plot Fx (s) in both PBE and PBEsol by dotted dashedand dotted curves.

1.0
1.2
1.4
1.6
1.8
2.0
2.2
2.6
2.7
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0

ccPBE
@Z = 1 Error(%)

ccPBE@Z = 2Error(%)

0.89459
0.676700.55895
0.489480.445270.417850.40152
0.38202
0.47532
0.44295
0.385590.36097
0.349170.34335
0.34100
0.34064
0.34105
0.34168
0.342210.885450.667290.54989
0.481770.440130.42272
0.41786
0.402670.48558
0.455740.400130.37486
0.36299
0.356820.35406
0.35361
0.35394
0.35467
0.35525

0.07
0.40
0.53
0.72
–1.05
–0.04
–0.12
0.01
–0.14
0.22
–0.62
0.55
0.05
–1.05
–1.16
–0.69
–0.28
0.20
0.95

–0.96
–1.00
–1.10
–0.87
–2.19
1.13
3.94
5.41
2.01
3.11
3.13
4.42
4.01
2.83
2.63
3.10
3.49
4.01
4.79

structure appears by reconciliation of the two constraints,
one from the gradient expansion, equation (29) or (34),and the other from the cusp condition which is explicitly
shown by “the cusp lines” in Figure 4.
Fig. 4. ∂Fx /∂s2 corresponding to Fx (s) in ccPBE in Figure 3.
It coincides with the target function (or the cusp line) at eachZ, plotted by the double-dotted-dashed curve, for s in therange 0.32−0.473.

parameter set so obtained is given in Table 2 in whichanother set of parameters fulfilling the cusp condition for Z = 2 and 3 are also added. Note that theparameters depending on Z are only those concernedwith F2 (s).
With those parameters, we can concretely give Fx (s)in ccPBE and its derivative ∂Fx (s)/∂s2 , both of which
are plotted in Figures 3 and 4, respectively. We haveplotted Fx (s) for three cases of the target functions
with Z = 1, 2, and 3, but its Z-dependence is found tobe weak. Compared with Fx (s) in PBE and PBEsol,
Fx (s) in ccPBE is enhanced much and has a characteristic structure for s < 0.6 but it increases smoothly fors > 0.6 and its actual value comes to the middle of PBE
and PBEsol.
As for ∂Fx (s)/∂s2 , its Z-dependence is much stronger
than that for Fx (s). The first sharp peak at s ≈ 0.2 isfound to be important to control the actual values of n(0)for rs < 2. (See the dependence of the peak structure on
Z in Fig. 4 and the change of n(0) with Z in Tab. 3.) This

3.3 On-top density in ccPBEIn ccPBE, we have successfully applied to the protonembedded electron-gas sphere and obtained very goodresults for n(0) in the entire range 1−12, as given
in Table 3. The TF-Kondo transition is seen by thejump in n(0) at rs ≈ 2.6. Accuracy of the results in
ccPBE is estimated by the relative error with respect
to the DMC data, given in %. Note that nineteen independent data in DMC are reproduced very well by
the appropriate choice of only five free parameters in
ccPBE.
For the target function with Z = 1, which agrees withthe atomic number of proton, the errors are at most about
1% but mostly much less than 1%. If we employ thetarget function with Z = 2, which is twice as large as
the atomic number of proton, in determining F2 (s) inequation (25) or (30), the errors are about several %,which may be said to be much larger than the case of
Z = 1, but at the same time it may be said to be still
small enough compared to the case of the original PBE.
This better performance may be said to be due to themuch smaller error in δZ/Z in ccPBE, even though we do
not employ the target function with the correct value of Z.
Incidentally, if we calculate n(0) in ccPBE with using thetarget function with Z = 3, the relative errors are foundPage 8 of 9

to be still not large, ranging from –1.3% to 7.0%, about
twice as large as those in the case of Z = 2.

4 Summary and discussion
By imposing the constraint originating from the cusp
theorem on the PBE scheme in GGA to DFT, wehave successfully constructed a new exchange-correlationenergy functional, referred to as ccPBE (cusp-corrected
PBE), and accurately reproduced the DMC data on the
on-top electron density n(0) in the proton-embedded electron gas with the density parameter rs in the range
1−12.
Five comments are in order:(i) Among fifteen parameters in the definition of Fx (s)in equation (25), only five parameters, namely, s0 ,s1 , s2 , µ1 , and B0 , can be chosen freely and independently of various constraints. After a rather extensive
search for appropriate values for them, we come tonotice that the adequate ranges for s1 and s2 arelimited by the cusp region (0.32, 0.473) in s-variable
space and probably the best values for them arethose in Table 2. In this sense, Fx (s) for s in therange s < 0.473 is almost completely determined
nonempirically by both the exact gradient expansion and the cusp theorem. As for other parameters,
namely, s0 , µ1 , and B0 having strong influence on
Fx (s) for s > 0.473, it is still not certain whether
the set of those values in Table 2 are best or not. A
better set of those parameters might be found in the
future.
(ii) As related to the above point, it might be considered that ρ0 in equation (33) is another independentand important parameter, but it does not seem to be
the case, because even if ρ0 = 1.0 is chosen insteadof ρ0 = 0.2 and consequently much different values
for Ci are used to define F2 (s), the self-consistently
determined results for n(r) do not change much,
indicating that we may choose any value of ρ0 as longas it is in the physically appropriate range 0.2−1.0.
(iii) As for the choice of Z in determining the target function, it is perfectly reasonable to choose Z = 1 forthe problems on hydrogen and the parameter set at
Z = 1 can be applied as it is to the solid hydrogen under high pressures. Even for the case of othervalues of Z, we might say that ccPBE with theparameter set at Z = 1 may provide better results
than PBE, but this needs to be confirmed by a comprehensive test of ccPBE for a wide class of realmaterials in the future. This test will also contributemuch to the choice of best appropriate values for theparameters s0 , µ1 , and B0 .
(iv) From a fundamental point of view, the xc functional should be universal and must be determined
only by the electron density n(r) itself. Thus one
may argue that the Z-dependent xc functional cannot be acceptable from the basic principles of DFT.
In order to overcome this criticism, we may propose the following amendment: Among three terms
in equation (25), only F2 (s) depends on Z throughEur. Phys. J. B (2018) 91: 189

the Z-dependence in C0 , . . ., C6 , and s2 . Then, let
us rewrite equation (30) as" i #
6X
p2
p
F2 (s, Z) =C0 (Z) +Ci (Z) 24s2 (Z)s2
i=1
× exp(−p/s2 (Z)2 ).

(38)

Now, since the term F2 (s, Z) becomes important
only in the cusp region at which the relation of Z =(9π/4)1/3 s/rs holds, we use its relation to introduce
the rs -dependent functional Fx (s, rs ) as
Fx (s, rs ) = F0 (s)+F1 (s)+F2 (s, (9π/4)1/3 s/rs ), (39)instead of Fx (s) in equation (25). This functional
Fx (s, rs ) satisfies the basic principles of DFT andat the same time the cusp theorem will be satisfied
for any Z. Note that once we consider the rs dependence in Fx , there are additional terms in vσx (r)in equation (15) and due changes must be made in
the subsequent calculations, including the determination of the coefficients Ci and their Z-dependence.
All those tasks concerning this amendment must be
done before implementing a comprehensive test of
ccPBE. Those works are left for the future.
(v) It is argued that the cusp theorem is satisfied inmeta-GGA [35,39]. Then the DMC data in Table 1
provide a good testing ground for meta-GGA. In
particular, it would be interesting to see which isthe predominant scheme among several proposed
ones [40–42,53–59] in reference to the DMC data.

Author contribution statement
YT as a single author has done all the works needed forthe preparation of this manuscript.

References
1. Z.D. Popovic, M.J. Stott, Phys. Rev. Lett. 33, 1164 (1974)
2. C.O. Almbladh, U. von Barth, Z.D. Popovic, M.J. Stott,Phys. Rev. B 14, 2250 (1976)
3. E. Zaremba, L.M. Sander, H.B. Shore, J.H. Rose, J. Phys.
F: Met. Phys. 7, 1763 (1977)
4. P. Jena, K.S. Singwi, R.M. Nieminen, Phys. Rev. B 17,
301 (1978)
5. G.W. Bryant, G.D. Mahan, Phys. Rev. B 17, 1744 (1978)
6. G.W. Bryant, Phys. Rev. B 19, 2864 (1979)
7. P. Jena, F.Y. Fradin, D.E. Ellis, Phys. Rev. B 20, 3543
(1979)
8. J.K. Nørskov, Phys. Rev. B 20, 446 (1979)
9. M.J. Stott, E. Zaremba, Phys. Rev. B 22, 1564 (1980)
10. M.J. Puska, R.M. Nieminen, M. Manninen, Phys. Rev. B
24, 3037 (1981)
11. M.J. Puska, R.M. Nieminen, Phys. Rev. B 27, 6121 (1983)
12. M.J. Puska, R.M. Nieminen, Phys. Rev. B 43, 12221(1991)
13. J.-H. Song, Ph. D. Thesis at Oregon State University, 2004,
http://hdl.handle.net/1957/29170

Eur. Phys. J. B (2018) 91: 189
14. V.U. Nazarov, C.S. Kim, Y. Takada, Phys. Rev. B 72,
233205 (2005)
15. Z.D. Popović, M.J. Stott, J.P. Carbotte, G.R. Piercy,
Phys. Rev. B 13, 590 (1976)
16. P. Jena, K.S. Singwi, Phys. Rev. B 17, 3518 (1978)
17. A.K. Gupta, P. Jena, K.S. Singwi, Phys. Rev. B 18, 2712
(1978)
18. P. Jena, A.K. Gupta, K.S. Singwi, Phys. Rev. B 18, 2723
(1978)
19. J. Arponen, E. Pajanne, J. Phys. C: Solid State Phys. 12,3013 (1979)
20. J. Gondzik, H. Stachowiak, J. Phys. C: Solid State Phys.
18, 5399 (1985)
21. L.M. Scarfone, A. Enver, Phys. Rev. B 43, 2272 (1991)
22. G. Sugiyama, L. Terray, B.J. Alder, J. Stat. Phys. 52, 1221
(1988)
23. A.I. Duff, J.F. Annett, Phys. Rev. B 76, 115113 (2007)
24. L.H. Thomas, Proc. Cambridge Philos. Soc. 23, 542 (1927)
25. E. Fermi, Rend. Accad. Naz. Lincei 6, 602 (1927)
26. P. Debye, E. Hückel, Phys. Z. 24, 185 (1923)
27. Y. Takada, R. Maezono, K. Yoshizawa, Phys. Rev. B 92,
155140 (2015)28. A.C. Hewson, Cambridge Studies in Magnetism, in The
Kondo Problem to Heavy Fermions, edited by D. Edwards,
D. Melville (Cambridge University Press, Cambridge,
England, 1993)
29. C.A. Kukkonen, A.W. Overhauser, Phys. Rev. B 20, 550(1979)
30. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77,
3865 (1996)[1pc]
31. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 78,
1396(E) (1997)
32. P.R. Antoniewicz, L. Kleinman, Phys. Rev. B 31, 6779
(1985)
33. J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov,
G.E. Scuseria, L.A. Constantin, X. Zhou, K. Burke, Phys.
Rev. Lett. 100, 136406 (2008)
34. J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov,
G.E. Scuseria, L.A. Constantin, X. Zhou, K. Burke, Phys.
Rev. Lett. 102, 039902(E) (2009)

Page 9 of 935. V.N. Staroverov, G.E. Scuseria, J. Tao, J.P. Perdew, Phys.
Rev. B 69, 075102 (2004)
36. T. Kato, Commun. Pure Appl. Math. 10, 151 (1957)
37. A.E. Carlsson, N.W. Ashcroft, Phys. Rev. B 25, 3474
(1982)
38. X.-Y. Pan, V. Sahnit, Phys. Rev. A 67, 012501 (2003)
39. J.P. Perdew, A. Ruzsinszky, J. Tao, V.N. Staroverov, G.E.
Scuseria, G. Csonka, J. Chem. Phys. 123, 062201 (2005)
40. J. Tao, J.P. Perdew, V.N. Staroverov, G.E. Scuseria, Phys.
Rev. Lett. 91, 146401 (2003)
41. J.P. Perdew, A. Ruzsinszky, G.I. Csonka, L.A. Constantin,
J. Sun, Phys. Rev. Lett. 103, 026403 (2009)42. J. Sun, A. Ruzsinszky, J.P. Perdew, Phys. Rev. Lett. 115,
036402 (2015)43. S. Azadi, W.M.C. Foulkes, Phys. Rev. B 88, 014115 (2013)
44. R.C. Clay III, M. Holzmann, D.M. Ceperley, M.A.
Morales, Phys. Rev. B 93, 035121 (2016)45. J.P. Perdew, Y. Wang, Phys. Rev. B 45, 13244 (1992)
46. S.-K. Ma, K.A. Brueckner, Phys. Rev. 165, 18 (1968)
47. Z. Wu, R.E. Cohen, Phys. Rev. B 73, 235116 (2006)
48. Z. Wu, R.E. Cohen, Phys. Rev. B 78, 197102 (2008)
49. Y. Zhao, D.G. Truhlar, Phys. Rev. B 78, 197101 (2008)
50. C.D. Hu, D.C. Langreth, Phys. Rev. B 33, 943 (1986)
51. P.S. Svendsen, U. von Barth, Phys. Rev. B 54, 17402(1996)
52. E.H. Lieb, S. Oxford, Int. J. Quantum Chem. 19, 427(1981)
53. Y. Zhao, D.G. Truhlar, J. Chem. Phys. 125, 194101(2006)
54. Y. Zhao, D.G. Truhlar, J. Phys. Chem. A 110, 13126(2006)
55. Y. Zhao, D.G. Truhlar, Theor. Chem. Acc. 120, 215(2008)
56. R. Armiento, S. Kümmel, Phys. Rev. Lett. 111, 036402(2013)
57. T. Aschebrock, R. Armiento, S. Kümmel, Phys. Rev. B
96, 075140 (2017)58. J. Tao, Y. Mo, Phys. Rev. Lett. 117, 073001 (2016)
59. Y. Mo, R. Car, V.N. Staroverov, G.E. Scuseria, J. Tao,Phys. Rev. B 95, 035118 (2017)