A simple correlation energy functional for the uniform electron gas is derived based on the second-order Moller-Plesset perturbation theory. It can reproduce the known correlation functional in the high-density limit, while in the mid-density range maintaining a good agreement with the near-exact correlation energy of the uniform electron gas to within 2 × 10−3 hartree. The correlation energy is a function of a density parameter rs and is of the form a * ln ( 1 + b r s + b r s 2 ) . The constants “a” and “b” are derived from the known correlation functional in the high-density limit. Comparisons to the Ceperley-Alder’s near-exact Quantum Monte Carlo results and the Vosko-Wilk-Nusair correlation functional are also reported.

Electron correlation energy functional is an important part of the density functional theory.1 In the local-density approximation of a uniform electron gas, general expression of the correlation energy is not known,2 except for the high- and low-density limits. There are many parameterization schemes such as the Vosko-Wilk-Nusair3 and the Perdew-Wang4 forms. These are primarily based on fitting functions to the near-exact Quantum Monte Carlo results of Ceperley and Alder.5 However, some of these functional forms such as the Vosko-Wilk-Nusair form are usually analytically complicated and inherently contain numerical irregularities.4 

In this report, a simple analytical function describing electron correlation energy within the local-density approximation is derived based on the second-order Moller-Plesset perturbation theory (MP2).6–8 The functional is of the form

ϵ c r s = a ln ( 1 + b r s + b r s 2 ) .
(1)

Not only it approaches the leading term9,10

ϵ c H D r s = A ln r s + C
(2)

in the known high-density limit but is also accurate in describing the correlation energy in the mid-density range. Here, r s = 4 π ρ / 3 1 / 3 is the Wigner-Seitz density parameter; ρ is the electron density. Atomic units are used throughout.

The constants “a” and “b” are

a = ln 2 1 2 π 2 , b = 20 . 456 255 7 .
(3)

As shown in Fig. 1, when the expression is plotted against the near-exact Quantum Monte Carlo results,3,5 it is able to describe the correlation energy in the mid-density range (2 ≤ rs ≤ 100) quite accurately, having the root-mean-squared error of only 1.856 × 10−3 hartree.

FIG. 1.

Comparisons to Monte Carlo results, VWN functional, and high-density limit.

FIG. 1.

Comparisons to Monte Carlo results, VWN functional, and high-density limit.

Close modal

The VWN functional also agrees with the Monte Carlo data as well. This is to be expected because the VWN parameters have been purposely fitted to these data points, whereas the correlation energy from Eq. (1) is derived independently. Interestingly, the root-mean-squared error of the VWN functional is 1.983 × 10−3 hartree, higher than that of Eq. (1). The detailed comparisons are provided in the supplementary material.11 

The form of the functional in Eq. (1) can be deduced as follows. From the expression of the MP2 correlation energy for a uniform electron gas6 

ε c 3 π 2 β 1 d q q 4 × 0 q 0 q d u d v u v 3 π 5 ρ 1 / 3 u + v 3 π 5 ρ 1 / 3 u ln u 2 + v ln v 2 ,
(4)

the integrand can be expanded as a polynomial in ρ1/3 to obtain the correlation energy at low density. For example,

ε c 3 π 2 β 1 d q q 4 0 q 0 q d u d v η 1 ρ 1 / 3 + η 2 ρ 2 / 3 + ,
(5)

where η1, η2 depend on the integrating variables u, v. Then, the integral can be separated into multiple terms, representing each order of ρ1/3 in the polynomial:

ε c 3 π 2 β 1 d q q 4 0 q 0 q d u d v η 1 u , v ρ 1 / 3 + 3 π 2 β 1 d q q 4 0 q 0 q d u d v η 2 u , v ρ 2 / 3 = η 1 ( I ) ρ 1 / 3 + η 2 ( I ) ρ 2 / 3 ,
(6)

where η 1 ( I ) , η 2 ( I ) are now constants.

However, Eq. (6) is not valid at high density because the expression does not approach infinity logarithmically as ρ → ∞. In other words, it is inconsistent with Eq. (2). The problem can be remedied by rewriting Eq. (6) inside a logarithmic function. Namely,

ε c = a ln 1 + η 1 ( I ) a ρ 1 / 3 + η 2 ( I ) a ρ 2 / 3 .
(7)

In the low-density range ρ1/3 ≪ 1, Eq. (7) falls back to Eq. (6) because ln 1 + x x if x is appreciably small. Furthermore, in the high-density range ρ1/3 ≫ 1, the term η 2 ( I ) a ρ 2 / 3 dominates; and ε c a ln η 2 ( I ) a ρ 2 / 3 = 2 a ln ρ 1 / 3 + Constant , which is essentially the same expression as Eq. (2). Therefore, Eq. (7) is also in its correct form in the high-density limit.

Note that, the density parameter rs is proportional to ρ−1/3. Therefore, Eq. (7) can be written in terms of rs as follows:

ϵ c = a ln 1 + b r s + c r s 2 ,
(8)

where “b” and “c” are constants. Equation (8) is valid when the two extreme cases are considered. When rs ≫ 1, the term 1 + b r s dominates, and ϵc is correct up to the leading term a b r s . On the other hand, when rs ≪ 1, the term c r s 2 is “turned on” instead, causing ϵc to approach the correct high-density original form. However, what happens in the mid-density range requires an assumption about the turning point at which electron gas “chooses” to behave like the low- or the high-density limit.

In the atomic unit, the Wigner-Seitz density parameter rs is defined as 4 π ρ / 3 1 / 3 . It defines the radius of a sphere which contains exactly one electron. In other words, it tells how far apart the electrons are from their nearest neighbors.

In the SI unit, the parameter rs is divided by the Bohr radius a0, the radius of a hydrogen atom, which can be viewed as the “comfortable” radius of an electron cloud encapsulating a unit of background positive charge. Therefore, the situation rs < 1 can be interpreted as the case when the density of electron gas is “squeezed” to be denser than that of a hydrogen atom. This causes the electrons to overlap and the Pauli exclusion principle to be significant, giving rise to the exchange interaction along with other quantum mechanical properties of the system.

When rs ≫ 1, the electrons are far apart. Each electron encapsulates one unit of background positive charge. Following Wigner’s view,12 there is plenty of room for an electron, and there is no need to overlap with its neighbor. In this limit, the system behaves classically, and the Coulomb interaction becomes the dominating term. It is no surprise that the correlation energy (per electron) is proportional to 1/rs in the case of the low-density limit.

Inspecting Eq. (8), one realizes that the parameters “b” and “c” govern the relative strength between the terms b r s and c r s 2 . For example, if b r s is the larger one, then the functional ϵc behaves like the low-density limit, and vice versa. At the tipping point between the low- and high-density limits, the two terms must be equal.

If it is assumed that such a tipping point occurs at rs = 1, exactly at the Bohr radius, then it follows that

b = c ,
(9)

as indicated in Eq. (1). The constants “a” and “b” in Eq. (1) are calculated by requiring that Eq. (1) approaches Eq. (2) as rs → 0. Namely, when rs is small, the term b r s 2 becomes the most dominant. In this case,

ϵ c r s a ln b r s 2 = 2 a ln r s + a ln b,
(10)

yielding a = A 2 , b = exp 2 C A . The values A and C have been analytically calculated10 using the perturbation theory. For example, in the paramagnetic case A = 1 ln 2 π 2 , C = 0 . 046 920 3 , giving the evaluated “a” and “b” in Eq. (3).

For a spin-polarized uniform electron gas, the correlation energy functional can be written as2 

ϵ c r s , ζ = ϵ c 0 r s + ϵ c 1 r s ϵ c 0 r s f ( ζ ) ,
(11)

where ζ ρ α ρ β ρ , ϵ c 0 r s , ϵ c 1 r s are the spin-polarization parameter, the paramagnetic-, and the ferromagnetic-correlation energy functional, respectively. The function f ζ = 1 + ζ 4 / 3 + 1 ζ 4 / 3 2 2 2 1 / 3 1 is a weighting factor between the two extreme cases: the paramagnetic and the ferromagnetic systems.

ϵ c 0 r s is essentially the same as Eq. (1). For the ferromagnetic ϵ c 1 r s , the values A and C from the high-density limit are needed to compute “a” and “b.” Using A = 1 ln 2 2 π 2 , C = 0 . 025 737 5 10 yields

a = ln 2 1 4 π 2 , b = 27 . 420 360 9 .
(12)

The correlation energy ϵ c 1 r s in the form given by Eq. (1) with the constants “a” and “b” listed above also agrees extremely well with the Monte Carlo results3,5 for the ferromagnetic case. The error in this case is well below 0.5 × 10−3 hartree as indicated in Fig. 3 of the supplementary material.11 Interestingly, this error is also comparable to the much more complicated formula obtained from an interpolation between high- and low-density limits.13 

In conclusion, a correlation energy functional is derived based on the second-order Moller-Plesset perturbation theory. The second-order energy correction is first expanded in powers of ρ1/3 to compute the correlation energy in the low-density range. Without losing significant accuracy, it is argued that when the expansion (up to the second order) is embedded inside a logarithmic function ln(1 + x), the analytical form of the high-density limit is recovered. Based on the assumption that the tipping point between the low- and high-density limits occurs when the Wigner–Seitz radius rs is exactly one Bohr radius, a simple analytic function ϵ c r s = a ln ( 1 + b r s + b r s 2 ) appears naturally. The constants “a” and “b” are then evaluated using the known high-density limit9,10 from the perturbation theory. The correlation energy functional agrees very well with the Monte Carlo results3,5 in the mid- and low-density range (2 ≤ rs ≤ 100), while approaching the known theoretical functional Alnrs + C in the high-density limit.

Its simple and accurate analytical form is useful as a starting point to further take into account the non-uniform density in molecules or bulk solids within the framework of the generalized gradient approximation (GGA).14–16 

1.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
2.
R. G.
Parr
and
W.
Yang
,
Density Functional Theory of Atoms and Molecules
(
Oxford University Press
,
Oxford
,
1989
).
3.
S. H.
Vosko
,
L.
Wilk
, and
M.
Nusair
,
Can. J. Phys.
58
,
1200
(
1980
).
4.
J. P.
Perdew
and
Y.
Wang
,
Phys. Rev. B
45
,
13244
(
1992
).
5.
D. M.
Ceperley
and
B. J.
Alder
,
Phys. Rev. Lett.
45
,
566
(
1980
).
6.
P. F.
Loos
and
P. M. W.
Gill
,
Int. J. Quantum Chem.
112
,
1712
(
2012
).
7.
C.
Møller
and
M.
Plesset
,
Phys. Rev.
46
,
618
(
1934
).
8.
G. S.
Handler
,
Int. J. Quantum Chem.
33
,
173
(
1988
).
9.
M.
Gell-Mann
and
K. A.
Brueckner
,
Phys. Rev.
106
,
364
(
1957
).
10.
P. F.
Loos
and
P. M. W.
Gill
,
Phys. Rev. B
84
,
033103
(
2011
).
11.
See supplementary material at http://dx.doi.org/10.1063/1.4958669 for additional figures and detailed comparisons to the VWN functional and the Quantum Monte Carlo data.
12.
E. P.
Wigner
,
Phys. Rev.
46
,
1002
(
1934
).
13.
J.
Sun
,
J. P.
Perdew
, and
M.
Seidl
,
Phys. Rev. B
81
,
085123
(
2010
).
14.
J. P.
Perdew
,
J. A.
Chevary
,
S. H.
Vosko
,
K. A.
Jackson
,
M. R.
Pederson
,
D. J.
Singh
, and
C.
Fiolhais
,
Phys. Rev. B
46
,
6671
(
1992
).
15.
A. D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
).
16.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).

Supplementary Material