Full PDF - Quantum Theory Project

PHYSICAL REVIEW A 90, 022504 (2014)
Comparative studies of density-functional approximations for light atoms in strong magnetic fields
Wuming Zhu
Department of Physics, Hangzhou Normal University, 16 Xuelin Street, Hangzhou, Zhejiang 310036, China
Liang Zhang
School of Information Science and Engineering, Hangzhou Normal University, 16 Xuelin Street, Hangzhou, Zhejiang 310036, China
S. B. Trickey
Quantum Theory Project, Department of Physics and Department of Chemistry, P.O. Box 118435, University of Florida,
Gainesville, Florida 32611-8435, USA
(Received 30 May 2014; published 11 August 2014)
For a wide range of magnetic fields, 0 B 2000 a.u., we present a systematic comparative study of
the performance of different types of density-functional approximations in light atoms (2 Z 6). Local,
generalized-gradient approximation (GGA; semilocal), and meta-GGA ground-state exchange-correlation (xc)
functionals are compared on an equal footing with exact-exchange, Hartree-Fock (HF), and current-densityfunctional-theory (CDFT) approximations. Comparison also is made with published quantum Monte Carlo
data. Though all approximations give qualitatively reasonable results, the exchange energies from local and
GGA functionals are too negative for large B. Results from the Perdew-Burke-Ernzerhof ground-state GGA
and Tao-Perdew-Staroverov-Scuseria (TPSS) ground-state meta-GGA functionals are very close. Because of
confinement, self-interaction error in such functionals is more severe at large B than at B = 0, hence selfinteraction correction is crucial. Exact exchange combined with the TPSS correlation functional results in a
self-interaction-free (xc) functional, from which we obtain atomic energies of comparable accuracy to those
from correlated wave-function methods. Specifically for the B and C atoms, we provide beyond-HF energies in
a wide range of B fields. Fully self-consistent CDFT calculations were done with the Vignale-Rasolt-Geldart
(VRG) functional in conjunction with the PW92 xc functional. Current effects turn out to be small, and the
vorticity variable in the VRG functional diverges in some low-density regions. This part of the study suggests
that nonlocal, self-interaction-free functionals may be better than local approximations as a starting point for
CDFT functional construction and that some basic variable other than the vorticity could be helpful in making
CDFT calculations practical.
DOI: 10.1103/PhysRevA.90.022504
PACS number(s): 31.15.E−
I. INTRODUCTION
Density-functional theory (DFT) [1] has been enormously
successful in electronic structure calculations of atoms,
molecules, surfaces, and solids. The success was fostered by
the availability of relatively accurate yet simple approximations to the exchange-correlation (xc) energy functional (and
the corresponding potential in the Kohn-Sham equation). Examples are the local density and local spin-density approximations (LDA, LSDA) and generalized-gradient approximations
(GGAs). For atoms and molecules, the literature of numerical
studies comparing the accuracy of such approximations to
that of wave-function or many-body-perturbation methods (socalled ab initio methods) is too large to cite with any fairness.
For atoms subjected to intense magnetic fields (B > 1
hartree a.u. = 2.3505 × 105 T), however, comparatively little
is known about the behavior and relative quality of various
DFT approximations [2–5]. What is known for physical atoms
is limited (as far as we can ascertain) to the LDA. Considerably
more investigation has been done on quantum dots, but, as
systems confined to two dimensions, they are not of direct relevance. Some time ago, two of us showed that various approximate functionals in DFT and current DFT (CDFT) [6] generate
substantial discrepancies compared to exact results for the
Hooke’s atom in an external field [7]. Since the density functional is, in principle, universal, those discrepancies are not
1050-2947/2014/90(2)/022504(14)
encouraging, hence a critical issue is the extent to which such
discrepancies occur in real three-dimensional (3D) atoms. That
issue was sharpened by the subsequent discovery of peculiar
v-representability properties for the two-electron system [8,9].
The present paper addresses the issue directly by systematic
calculation of the energetics of light atoms (2 Z 6, Z =
atomic number) over a very large range of B (0 B 2000 a.u.) with extremely rich, aspherical basis sets and several
xc functionals, including one combination which apparently
has gone unexplored heretofore. We take advantage of the
availability of several Hartree-Fock (HF) [10–18] and correlated wave-function studies of at least a few light atoms in an
external B field [19–27] and of quantum Monte Carlo (QMC)
studies [28,29] for comparisons. In particular, comparison of
our HF calculations with prior ones provides calibration of the
accuracy of our techniques, while comparison with the QMC
results gives a measure of absolute accuracy and trends in
accuracy for the various DFT approximations.
Though the manifestly gauge-invariant current CDFT was
established more than two decades ago [6], CDFT calculations
on 3D systems remain rare, with the majority being for
vanishing or low-B-field cases [30–33]. There has been recent
work on molecular magnetic properties in CDFT up to about
B = 1 a.u. by Tellgren et al. [34]. In several ways their work
is complementary to ours, thus some of their findings and
technical details are discussed below.
022504-1
©2014 American Physical Society
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
CDFT parametrizes the many-electron ground state via the
electron number density n(r) and paramagnetic current density
jp (r). The current-independent limit is ordinary DFT with
appropriate explicit dependence on B, which is the second
form of DFT studied here. Because it is a simple extension of
existing zero-field DFT approximations to B = 0, that second
form enables better comprehension of the challenges in CDFT
calculations.
As suggested above, some of the approximate xc functionals are well tested in the B = 0 case [35–37], but their
performance for B = 0 is not. We did not include clusters and
solids for several reasons: High-level correlated results for
such systems are rare, there is the technical burden of developing accurate, fast multicenter integrals for the special basis
functions suitable for high B fields, and the study of atoms
avoids issues of gauge invariance and periodicity breaking.
The essentials of the DFT and CDFT formulations are in the
next section. Basis sets and numerical methods are addressed
in Sec. III (and the Appendix), followed by results in Sec. IV,
and a concluding summary in Sec. V.
The xc contributions are
+
1
1
.
2 i=j |ri − rj |
(1)
Here ri , mi , and ms,i are the space coordinate, magnetic
quantum number, and spin z component for the ith electron,
respectively. (As indicated already, Hartree atomic units
= melectron = qelectron = 1 are used throughout unless stated
explicitly; one Hartree a.u. of B field corresponds to 2.3505 ×
105 T.) Without explicit exception, unpaired electrons always
are taken as spin down in what follows, i.e., ms,i = − 12 for
unpaired electrons.
In CDFT, the variational minimization for the manyelectron ground state is mapped to a Kohn-Sham (KS) system
which generates n(r) and jp (r). The system xc energy thus
depends on both in general. The CDFT KS equation is
B2 2
1
Z
B
(x + y 2 ) + (mi + 2ms,i ) + vH (r)
− ∇2 − +
2
r
8
2
1
+ vxc (r) + Axc (r) · ∇ φi (r) = i φi (r),
(2)
i
with
vH (r) =
n(r )
dr .
|r − r |
δExc [n(r),jp (r)]
δn(r)
(4)
Axc (r) =
δExc [n(r),jp (r)]
,
δjp (r)
(5)
and
while the densities are
n(r) =
|φi (r)|2
(6)
i
and
jp (r) =
1 ∗
[φi (r)∇φi (r) − φi (r)∇φi∗ (r)].
2i i
The system total energy is
Etot = Ts + J + Exc [n,jp ] +
dr n(r) −
(7)
Z
r
B
B2 2
2
(x + y ) +
(mi + 2ms,i ),
+
8
2 i
II. METHODOLOGY
Upon imposition of a uniform, static external magnetic
field along the z direction, B = B zˆ , a central field atom
at the origin becomes cylindrically symmetric. The system
Hamiltonian commutes with rotations about the field direction,
so the magnetic quantum number m is still good. Three other
conserved quantities are the total spin S2 , its z component Sz ,
and spatial parity in z, hence a state can be labeled by those
quantum numbers. In the Coulomb gauge, A(r) = 12 B × r,
with A the vector potential, the system Hamiltonian then is
1
B
Z
B2 2
2
2
xi + yi + (mi + 2ms,i )
H=
− ∇i − +
2
ri
8
2
i
vxc (r) =
(8)
where Ts denotes the noninteracting (KS) kinetic energy and
J is the Hartree (classical electron-electron repulsion) energy.
In the spirit of linear response, it often is assumed that the
Axc (r) contribution to Etot is small compared to the ordinary
DFT Exc . The zeroth-order approximation to Exc [n,jp ] thus
takes the same form as the xc functional in ordinary DFT,
Exc [n], and the KS equation reduces to the familiar form except
that explicit external field contributions appear. We refer to
this approximation as naive B-DFT (nB-DFT). Observe that
in the nB-DFT approximation, only the explicit dependence
of Exc [n,jp ] upon jp (r) is neglected, that is, Exc [n,jp ] ≈
Exc [n,0]. Other external field dependencies, including the
interaction between jp and the external B field, and the implicit
dependence on jp (r), are unaltered.
The diversity of approximate xc functionals for the B = 0
case is well documented. Following the Perdew-Schmidt
Jacob’s ladder of functional complexity [38], the first two
rungs, LDA and GGA, are local and semilocal functionals.
The third, meta-GGA (mGGA), has an explicit orbital dependence from the KS kinetic energy density. We choose
the PW92 [39], Perdew-Burke-Ernzerhof (PBE) [40], and
Tao-Perdew-Staroverov-Scuseria (TPSS) [41] functionals as
representatives of the respective rungs. As is well known, both
local and semilocal functionals suffer from improper electron
self-interaction. Among the correlation functionals for those
three, only the TPSS correlation functional is self-interaction
free [41,42]. Therefore, the combination of exact-exchange
(EXX) and the TPSS correlation functional,
Exc = ExXX + EcTPSS ,
(9)
is a self-interaction-free xc functional. Exact exchange in DFT
[43] is defined in terms of occupied KS orbitals {φi (r)} in the
Fock integral,
∗
φi (r)φj∗ (r )φj (r)φi (r )
1
XX
dr dr .
Ex [{φi (r)}] = −
2 i,j
|r − r |
(3)
022504-2
(10)
COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . .
PHYSICAL REVIEW A 90, 022504 (2014)
Since ExXX is fully nonlocal, the combination in Eq. (9) corresponds to the fourth rung of the Jacob’s ladder or hyper-GGA
(HGGA). Note that EcTPSS is not fully compatible with ExXX ,
since the TPSS correlation hole is rather local and misses the
long-range part. Because the spatial extension of atomic systems usually is limited such that long-range hole components
are negligible or nearly so [44], this combination should not
be a problem for the present study. That appraisal is confirmed
numerically (see below). However, one cannot assume that
Eq. (9) also would be generally applicable to extended systems,
such as molecules or solids. For a correlation functional fully
compatible with exact exchange, we would need to move
up to the fifth rung of the Jacob’s ladder. An example is
the generalized random phase approximation (RPA) [45].
However, such RPA calculations would involve significant
additional development, a task we leave to a future study.
The implicit dependence of Exc on jp in the nB-DFT approximation does not matter for local and semilocal functionals, but does for orbital-dependent functionals. Equation (10)
is a specific case. Since KS orbitals depend on jp , and ExXX
in turn depends on those KS orbitals, there is an indirect
dependence of ExXX on jp . Thus, the current contribution
to the exchange energy is included at least partially in the
nB-DFT approximation with exact exchange [see Eq. (9)]. An
additional complication is that to stay strictly within the KS
scheme therefore would require solution of the corresponding
optimized effective potential (OEP) equations [46]. The OEP
method and associated simplifying approximations have been
generalized to the current-spin-DFT (CSDFT) case [47].
However, the numerical burden of the OEP equations is
substantial, with the result that direct variation of Eq. (9)
with respect to the orbitals frequently is used instead. Called
“generalized KS” in the quantum chemistry literature, this
approach gives one-electron equations of the same structure as
the HF equations but with a local correlation potential added.
Numerical results show that atomic energies obtained from
OEP and generalized KS calculations are quite close for the
X-only functional [46,48]. Thus we expect the energies also
to be close for OEP and generalized KS calculations when
Eq. (9) is used, because its correlation part is a mGGA. For
computational simplicity, we therefore used the generalized
KS procedure for orbital-dependent functionals.
In contrast to the nB-DFT situation, there is little choice
in xc functionals for CDFT. To our knowledge, the only one
published for 3D systems that is in a form implementable
for self-consistent solution of the KS equations is the local
approximation due to Vignale, Rasolt, and Geldart (VRG)
[6]. From the perturbative energy of a homogeneous electron
gas (HEG) in uniform B and gauge-invariance arguments,
Vignale and Rasolt gave the local approximation for the current
contribution to the xc energy [6],
VRG
Exc
[n(r),ν(r)] =
g(r) |ν(r)|2 dr,
(11)
where
jp (r)
,
n(r)
χ (n(r))
kF
−
1
.
g(r) = g(n(r)) =
24π 2 χ0 (n(r))
ν(r) ≡ ∇ ×
(12)
(13)
The variable ν is the vorticity mentioned already, while χ and
χ0 are the orbital magnetic susceptibilities for the interacting
and noninteracting HEG, respectively. kF = (3π 2 n)1/3 is the
usual local Fermi momentum.
From computed susceptibility data over the range 1 rs 10 tabulated in Ref. [49], Lee, Colwell, and Handy (LCH)
generated an analytical fit for the ratio [30],
χ
,
(14)
s≡
χ0
sLCH = (1.0 + 0.028rs )e−0.042rs
(15)
3 1/3
with rs = ( 4πn
) the usual Wigner-Seitz radius. The resulting
approximation is
gLCH =
kF
(sLCH − 1) .
24π 2
(16)
With various regularizations (cutoffs) which we discuss
shortly, the VRG functional has been used to calculate
magnetizabilities [30,34], nuclear shielding constants [31], and
frequency-dependent polarizabilities [32] for small molecules,
and ionization energies for atoms [33]. Because all those
studies are for B → 0, they are inconclusive regarding the
general, B 0 utility of the VRG functional. Except for
our Hooke’s atom results, Ref. [7], little is known about the
field-dependent behavior of finite systems for very large B that
emerges from the VRG functional.
An important technical issue is evident in Eq. (15), namely,
the sensitivity of the predicted physics to the details of the
low-density extrapolation. Orestes, Marcasso, and Capelle
(OMC) proposed two other fits, both polynomial [33]. A fit
by Tao and Perdew [50] actually gives the correct rs → 0 limit
and physically plausible low-density behavior. Subsequently
a refinement was given by Tao and Vignale [51]. As discussed
in detail in Ref. [34], those two seem to differ little in
practice. That reference also shows how the LCH, Tao-Perdew,
and Tao-Vignale parametrizations tend to zero as rs → ∞,
whereas the OMC and VRG forms do not. They also studied
the consequences of the peculiar limiting behavior of the rather
different Higuchi and Higuchi [52] functional.
The relevant point here is that even at the comparatively
low fields considered in Ref. [34], there is considerable
variation in the numerical stability for that collection of
low-density regularizations but not large shifts in the computed
quantities. Thus we resorted to work we had done earlier on a
different regularization [7]. From a study of many numerical
experiments, we arrived at
gcutoff =
kF
(c1 + c2 rs )e−αcutoff rs
24π 2
(17)
for densities below ncutoff . The values of c1 and c2 follow
from requiring a smooth connection between gcutoff and g(n)
at the specified transition density ncutoff . In this work, we
use ncutoff = 0.001 bohr−3 and αcutoff = 2.0 bohrs−1 , unless
otherwise explicitly specified. The numerical value of ncutoff
corresponds to rs = 6.2 bohrs, so that the regularization
afforded by gcutoff could be tested and calibrated against the
available data on 6.2 rs 10 bohrs.
022504-3
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
where
III. BASIS SETS AND NUMERICAL METHODS
In molecular calculations, the KS equations usually are
solved by expansion in a Gaussian-type orbital (GTO) basis.
Less commonly, Slater-type orbitals are used, but they are
not well suited for systems in a strong B field. Specifically
Jones, Ortiz, and Ceperley used GTOs to calculate small
atoms and molecules in strong fields, and found that basis
elements with very large angular momentum (lmax up to 35)
are required for reasonable numerical convergence [15,53].
Tellgren et al. [34] used rich, conventionally isotropic GTO
basis sets and commented that their calculations thus were
limited to maximum B ≈ 1 a.u. For vanishing B field, Porezag
and Pederson introduced optimized GTO basis sets for DFT
calculations [54]. The difficulty for our purposes is that their
technique is based on total energy optimization and therefore
requires specification of the theoretical framework (HF or
DFT). Our intent is to compare the consequences of different
approximations and methods with the same basis sets for all
calculations in order to avoid artifacts.
Here we choose anisotropic-GTO (AGTO) basis sets such
as were introduced by Aldrich and Greene [55], and exploited
extensively by Schmelcher and Cederbaum [56]. In cylindrical
coordinates (ρ,z,φ), the basis functions are
χj (ρ,z,φ) = Nj ρ nρj znzj e−αj ρ
2
−βj z2 imj φ
e
,
where nρj = |mj | + 2kj , kj = 0,1, . . . , with mj = . . . ,
−2, −1,0,1,2, . . . , and nzj = πzj + 2lj ,lj = 0,1, . . . , with
πzj = 0,1.
AGTO basis sets have been used in several studies of
atoms in strong B fields [20,21,24,25,27] because of their
flexibility in describing elongation of the electron orbitals and
densities along the B-field direction. With no unique way
to determine the basis exponents, αj and βj , a physically
sound prescription must be developed. Kappes and Schmelcher
gave an optimization algorithm [57] which, however, must be
executed for each combination of atomic configuration and
field strength, a tedious, time-consuming task. Furthermore, a
good initial guess is required [21]. An alternative is the use of
nearly optimized basis sets. Kravchenko and Liberman (KL)
[58] investigated one-electron systems, the hydrogen atom and
the hydrogen molecular ion, and showed that systematically
constructed AGTO basis sets could provide accuracy of 10−6
hartree or better.
As detailed in the Appendix, we constructed KL-like
highly optimized basis sets without needing full nonlinear
optimization. Each is comprised of a sequence of subsets
related by
αj, = βj + (1 + μ )(B),
μ2,3 = ±0.2,
and, in practice, γ = B/Z 2 . The base sequence is μ=1 = 0.
Others are μ = ±0.2,±0.4 for the second, third, fourth,
and fifth sequences, respectively. For = 2,3, there are half
as many exponents (with doubled spacing) as in the base
sequence, while for = 4,5 there are one-fourth as many (with
quadrupled spacing) as in the base sequence. Similar to the KL
basis sets, we used even-tempered Gaussian (ETG) sequences
for the longitudinal exponents βj . Again, see the Appendix.
The kinetic energy, nuclear-attraction, overlap, and onecenter Coulomb repulsion matrix elements with respect to the
AGTO basis all can be expressed in closed forms, though the
Coulomb repulsion expression is lengthy [21]. The vorticity,
ν(r), was calculated analytically. For CDFT, both xc potentials,
vxc (r) and Axc (r), were evaluated numerically on a twodimensional mesh, with up to 500 points along each direction.
The convergence of total energies with respect to basis set size
and with respect to the number of mesh points was checked
carefully. The self-consistent field (SCF) convergence criterion
was 10−8 in total energy. In general the mesh convergence was
no worse than a few units in the last digit of the tabulated data.
j = 1,2,3, . . . ,
(18)
μ1 = 0.0,
b(γ ) = −0.16[tan−1 (γ )]2 + 0.77 tan−1 (γ ) + 0.74, (21)
μ4,5 = ±0.4,
(19)
with αj,1 ≈ βj ∀j and (B) a parameter optimized for each B.
The result of extensive numerical exploration (see Appendix)
is
B
4 βj −1/2
4 βj −2
(B) =
+ 1+
,
4 1+
20
b(γ ) B
b(γ ) B
(20)
IV. RESULTS AND DISCUSSION
A summary of our results at four B magnitudes (42.54 a.u. = 1.0 × 107 T, 212.72 a.u. = 5.0 × 107 T,
425.43 a.u. = 1.0 × 108 T, 2127.2 a.u. = 5.0 × 108 T) with
comparison to a broad selection of near-exact calculations,
is given in Tables I and II. For illustration, a selection of
detailed results is in Tables III and IV for the He and Li
atoms, again with comparison to published data. Those four
tables are the main focus of the discussion in this section.
More extensive tabulations are in the Supplemental Material
[59], which includes more electronic states for the He, Li, Be,
B, and C atoms and the positive ions, Li+ , Be+ , and B+ , in
field strengths 0 B 2000 a.u. In all cases, the electronic
states are labeled according to their corresponding zero-field
Hartree-Fock electronic configurations.
A. Unrestricted Hartree-Fock calculation with AGTO basis
We first did unrestricted Hartree-Fock (UHF) calculations
mainly as a check on the basis set construction and for
verification of numerical technique. The UHF energies are
in excellent agreement with published values. For the He
atom, Zhao and Stancil expanded the 1s 2 electron orbitals in
a B-spline basis and used quadruple-precision calculations to
obtain very accurate HF energies for 0 B 100 a.u. [17].
Our results differ from theirs by no more than a couple of
μ-hartree. For other electronic states, we compare with the
data by Jones, Ortiz, and Ceperley [20], who also used AGTO
basis sets. Frequently our energies are slightly lower, perhaps
because of better basis optimization.
Comparison data for atoms and ions with Z 3 are
from the studies by Ivanov and Schmelcher [16], who used
two-dimensional (2D) mesh methods. Differences usually are
less than 0.1 m-hartree, and the overall agreement is quite
satisfactory. Exceptions are the 1s 2 2s 2 state of the Be atom
022504-4
COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . .
PHYSICAL REVIEW A 90, 022504 (2014)
TABLE I. Atomic energies for He, Li, Be, B, and C atoms in B fields by different methods (energies in hartree, B in tesla). Present
work is “HF, pres.” plus the four xc functionals PW92 (LDA), PBE (GGA), TPSS (MGGA), and exact exchange plus TPSS correlation
(HGGA). Comparison data are for various QMC methods denoted as in Refs. [28] Tables I–IV for “RPDQMC” through “MCPH3” and [29]
for FPDQMC2.
Z B (T)
HF, pres.
LDA
GGA
MGGA
HGGA
RPDQMC FPDQMC VQMC
HFFEM
2DHF
MCPH3 FPDQMC2
2 1.0 × 107 −9.6969 −9.6304 −9.8332 −9.7988 −9.7216 −9.735 −9.735 −9.617 −9.393 −9.697 −9.603
5.0 × 107 −16.928 −17.151 −17.618 −17.535 −16.957 −17.00 −17.01 −16.88 −16.72 −16.93 −16.78
1.0 × 108 −21.314 −21.873 −22.544 −22.425 −21.344 −21.41 −21.41 −21.28 −21.15 −21.31 −21.19
5.0 × 108 −35.349 −37.790 −39.345 −39.074 −35.381 −35.51 −35.54 −35.35 −35.24a −35.35 −35.18
−9.714
−16.968
−21.374
−35.489
3 1.0 × 107
5.0 × 107
1.0 × 108
5. × 108
−19.61
−35.01
−44.61
−76.36
−19.891
−35.420
−45.107
−77.054
4 1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
−33.013 −32.916 −33.239 −33.126 −33.076 −33.18 −33.15 −32.50 −31.14 −33.01 −32.70
−59.395 −59.860 −60.571 −60.354 −59.471 −59.61 −59.57 −59.02 −58.14 −59.40 −58.80
−76.184 −77.301 −78.313 −78.036 −76.264 −76.47 −76.44 −75.89 −75.11 −76.18 −75.56
−132.68 −137.56 −139.86 −139.37 −132.76 −133.40 −133.25 −132.55 −131.97 −132.68 −131.78
−33.059
−59.504
−76.344
−133.067
5 1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
−48.961 −48.853 −49.240 −49.116 −49.044 −49.17 −49.17 −47.96 −45.53 −48.96 −48.77
−88.602 −89.189 −90.015 −89.765 −88.702 −88.97 −88.86 −87.94 −86.32 −88.60 −87.83
−114.26 −115.65 −116.82 −116.50 −114.37 −114.73 −114.62 −113.74 −112.34 −114.26 −113.37
−202.10 −208.14 −210.78 −210.25 −202.22 −203.04 −203.00 −201.79 −200.83 −202.10 −201.24
−49.021
−88.744
−114.467
−202.605
6 1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
−67.580 −67.467 −67.921 −67.796 −67.684 −67.88 −67.95 −65.85 −62.00 −67.58 −68.13
−122.61 −123.32 −124.26 −124.00 −122.73 −123.07 −122.93 −121.49 −118.92 −122.61 −121.57
−158.75 −160.42 −161.74 −161.40 −158.88 −159.42 −159.16 −157.80 −155.67 −158.75 −157.58
−284.26 −291.45 −294.41 −293.86 −284.41 −285.39 −285.32 −283.96 −282.20 −284.26 −282.78
−67.655
−122.783
−159.010
−284.874
a
−19.860
−35.345
−44.997
−76.781
−19.777
−35.691
−45.841
−80.471
−20.040
−36.283
−46.688
−82.417
−19.957
−36.122
−46.475
−82.010
−19.904
−35.397
−45.052
−76.840
−19.92
−35.50
−45.20
−77.28
−19.93
−35.47
−45.16
−77.21
−19.58
−35.19
−44.83
−76.73
−18.99
−34.76
−44.50
−76.44
−19.86
−35.35
−45.00
−76.78
An obvious sign error in the original is corrected here.
and B+ ion, and the 1s 2 2s 2 2p−1 state of the B atom, especially
in very large B fields [59]. Large discrepancies (from 0.5
mhartree at B = 1 a.u. to 1.9 hartree at B = 2000 a.u. for
the Be atom, for example) seem plausibly to be attributable
to different allowed spatial symmetries for the electron wave
function. They adopted asymmetrical wave functions for 2s 2
electrons with respect to the z = 0 plane [16], whereas we
required definite z parity, πz = 0,1. Excluding those three
instances, we can compare our HF energies with theirs
evenhandedly, because they did not find symmetry breaking
in any other cases. Usually our HF energies are slightly
higher than theirs, with differences no more than a few units
in the last quoted digit. Presumably those differences arise
from numerical mesh error in their calculations and basis set
truncation error in ours. On variational grounds, the basis set
truncation error always is positive, so our data are an upper
bound for the HF energies. In contrast, there can be small
magnitude but negative numerical errors from Ivanov and
Schmelcher’s 2D HF mesh method. Compare, for example,
their H atom results with those from Ref. [60]. Hence, one may
speculate that the true HF energies lie between our data and
theirs. The fact that our basis-based calculation is comparable
to the best numerical results in accuracy suggests strongly that
our construction should result in sufficiently large and flexible
basis sets for describing atoms in a wide range of B fields for
a range of theoretical methods.
B. Results of nB-DFT approximations
The PW92, PBE, and TPSS functionals (LDA, GGA,
and MGGA, respectively) were implemented and tested in
the context of the nB-DFT approximation discussed above.
Results again are listed in Tables I, III, and IV, and in the
Supplemental Material [59].
Regarding context, there are fewer published DFT calculations for atoms in B than for HF. Comparison of nonzero
B-field DFT results moreover is handicapped by the different
magnetic field values for which various authors present their
results, and by the different xc functionals implemented. The
functional due to Jones [2], which is at the level of LDA,
was used by Neuhauser, Koonin, and Langanke [11], and
by Braun [4]. The simple Dirac X-only functional was used
by Relovsky and Ruder [3], and by Braun [4]. Our results
for the Dirac X functional agree well with Braun’s. Medin
and Lai considered only the high-field case. Based on the
adiabatic approximation [5], they chose the Danz-Glasser
X functional and a strong-field-limit expression for the C
functional [61]. Both are local density approximations which
are somewhat crude in comparison to modern ground-state,
B = 0 xc functionals. And, it is not clear that a well-behaved
low-field limit of that xc functional exists or, if it does, that it
recovers the ordinary LDA. As far as we can tell, there is no
prior application of density-gradient-dependent functionals to
atoms in a strong B field. In the broader context, since DFT
calculations in principle include electron correlation, it also is
appropriate to compare our nB-DFT results with those from
correlated wave-function methods [19–27].
We address some of the more predictable outcomes first.
Inspection of Tables I, III, and IV shows that total atomic
energies from several different methods differ from a few hundredths of a hartree at low B up to a few hartrees at very large
B field. To illustrate the point, Table II recasts the results in
022504-5
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
TABLE II. Deviations from reference atomic energies (last column, FPDQMC2) for He, Li, Be, B, and C atoms in B fields by different
methods (energies in hartree, B in tesla). Present work is “HF, pres.” plus the four xc functionals PW92 (LDA), PBE (GGA), TPSS (MGGA),
and exact exchange plus TPSS correlation (HGGA). Comparison data are for various QMC methods denoted as in Refs. [28] Tables I–IV for
“RPDQMC” through “VQMC” and [29] for FPDQMC2.
Z
B (T)
HF, pres.
LDA
GGA
MGGA
HGGA
RPDQMC
FPDQMC
VQMC
2
1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
0.017
0.040
0.060
0.140
0.084
−0.183
−0.499
−2.301
−0.119
−0.650
−1.170
−3.856
−0.085
−0.567
−1.051
−3.585
−0.008
0.011
0.030
0.108
−0.021
−0.032
−0.036
−0.021
−0.021
−0.042
−0.036
−0.051
0.097
0.088
0.094
0.139
−9.714
−16.968
−21.374
−35.489
3
1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
0.031
0.075
0.110
0.273
0.114
−0.271
−0.734
−3.417
−0.149
−0.863
−1.581
−5.363
−0.066
−0.702
−1.368
−4.956
−0.013
0.023
0.055
0.214
−0.029
−0.080
−0.093
−0.226
−0.039
−0.050
−0.053
−0.156
0.311
0.230
0.277
0.324
−19.891
−35.420
−45.107
−77.054
4
1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
0.046
0.109
0.160
0.387
0.143
−0.356
−0.957
−4.493
−0.180
−1.067
−1.969
−6.793
−0.067
−0.850
−1.692
−6.303
−0.017
0.033
0.080
0.307
−0.121
−0.106
−0.126
−0.333
−0.091
−0.066
−0.096
−0.183
0.559
0.484
0.454
0.517
−33.059
−59.504
−76.344
−133.067
5
1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
0.060
0.142
0.207
0.505
0.168
−0.445
−1.183
−5.535
−0.219
−1.271
−2.353
−8.175
−0.095
−1.021
−2.033
−7.645
−0.023
0.042
0.097
0.385
−0.149
−0.226
−0.263
−0.435
−0.149
−0.116
−0.153
−0.395
1.061
0.804
0.727
0.815
−49.021
−88.744
−114.467
−202.605
6
1.0 × 107
5.0 × 107
1.0 × 108
5.0 × 108
0.075
0.173
0.260
0.614
0.188
−0.537
−1.410
−6.576
−0.266
−1.477
−2.730
−9.536
−0.141
−1.217
−2.390
−8.986
−0.029
0.053
0.130
0.464
−0.225
−0.287
−0.410
−0.516
−0.295
−0.147
−0.150
−0.446
1.805
1.293
1.210
0.914
−67.655
−122.783
−159.010
−284.874
Table I in the form of deviations of calculated values for various
approximations with respect to the FPDQMC2 values from
Ref. [29]. (There are not enough tabulated near-exact results
for a similar recasting of Tables III and IV, but study shows
similar results when comparison is possible.) Even within
DFT there is similar variation for different xc functionals.
All give similar, qualitatively correct results. The changes in
energies with increasing B field also are shown in Figs. 1 and 2
for the 1s2p−1 state of He and the 1s2p−1 3d−2 state of Li,
respectively. Clearly, the TPSS (meta-GGA) and PBE (GGA)
results are very close; the two curves are essentially indistinguishable. Moving from LDA to GGA significantly improves
the accuracy of atomic energies at low B field, but shifting from
GGA to MGGA (hence including the kinetic energy density
as an additional ingredient) does not offer a comparable gain.
Also see tables in the Supplemental Material [59].
Published results for correlated wave-function methods
(WFMs) are only for 0 B 100 a.u. for He and Li, and
for the Be atom, 0 B 10 a.u. From Figs. 1 and 2, it is
hard to tell which methods of those we are comparing agree
better with those data points (indicated by “+”). All curves
essentially are superimposed upon the quoted data.
It is easy to see why different methods do not give
drastically different outcomes. Equations (1) and (2) show that
the B field enters the system Hamiltonians in the same way
for all the methods. The field only modifies the single-electron
operators, with no effect on the two-electron operator. Effects
of a B field on the electron-electron interaction can occur only
indirectly through changing the many-electron wave function
or, in DFT, the KS orbitals. Since all methods handle the oneelectron part exactly and it contributes a large portion of the
total energy, results from different methods are not expected
FPDQMC2
to differ too much as long as the two-electron energy is treated
reasonably. Figure 3 illustrates these points by giving separated
one-electron and two-electron energies in the He atom triplet
state, 1s2p−1 . The curves for one-electron energies, labeled by
E1e , are nearly indistinguishable by method. Differences arise
mainly from the two-electron energies, E2e , precisely because
of differing evaluations of the Hartree-plus-xc energy, EHxc .
Detailed examination of the contributions to EHxc leads,
however, to what seem to be previously undocumented results.
Of course, both the Hartree term and the X term are much larger
(in magnitude) than the C energy. In HF calculations, electron
self-interaction energy ESI is canceled exactly by the selfexchange part of Ex and Ec is totally neglected. In DFT calculations, with either local or semilocal approximate xc functionals, the cancellation is incomplete, causing the so-called
self-interaction error. (We discuss only integer self-interaction
error, which is unexplored in the large B-field context, and
leave noninteger self-interaction error (SIE) [62] to separate
study.) Nevertheless, some gradient-dependent functionals do
behave remarkably well in some circumstances. An example
is the He atom low-field ground state 1s 2 . Define the selfinteraction energy in terms of the KS orbital self-repulsions as
1
2 i
N
ESI :=
|φi (r)φi (r )|2
dr dr .
|r − r |
(22)
Exact self-interaction cancellation for a two-electron singlet
is Ex [n]/J [n] = − 12 (or ESI /Ex = −1). The PW92 LDA X
functional achieves 87% cancellation for this He state at B =
0, compared to more than 99% for B < 0.4 a.u. for both the
PBE and TPSS X functionals. This approximate cancellation
breaks down completely when B > 10 a.u. (see Fig. 4). Both
022504-6
COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . .
PHYSICAL REVIEW A 90, 022504 (2014)
TABLE III. Atomic energies of the helium atom in B fields by different methods. (Energies in hartree, B in a.u. The xc functionals are
PW92 (LDA), PBE (GGA), and TPSS (mGGA), respectively. HGGA comprises exact exchange with the TPSS correlation functional.)
State
B (a.u.)
HF, present
0
0.08
0.1
0.5
0.8
1
2
5
8
10
20
50
80
100
800
1000
−2.861679
−2.860417
−2.859709
−2.814450
−2.746839
−2.688884
−2.289144
−0.532442
1.591275
3.110634
11.319611
38.143906
66.092087
85.004179
770.54396
968.44540
0
0.08
0.1
0.5
0.8
1
2
5
8
10
20
50
80
100
800
1000
−2.131347
−2.236463
−2.259234
−2.615549
−2.830207
−2.959686
−3.502049
−4.617248
−5.400409
−5.829510
−7.427702
−10.264491
−12.101321
−13.076652
−26.126547
−28.032093
0
0.08
0.1
0.5
0.8
1
5
8
10
50
80
100
800
1000
−2.055211
−2.166315
−2.187305
−2.500874
−2.687529
−2.800387
−4.276634
−4.987052
−5.378085
−9.455332
−11.154700
−12.058706
−24.229300
−26.015249
HF, lit.a
LDA
GGA
MGGA
HGGA
−2.83445
−2.83308
−2.83232
−2.78378
−2.71240
−2.65177
−2.23932
−0.45536
1.68598
3.21416
11.44959
38.29203
66.23035
85.13000
770.0003
967.7239
−2.89294
−2.89160
−2.89085
−2.84347
−2.77363
−2.71423
−2.30882
−0.54618
1.57640
3.09325
11.28092
38.01980
65.87939
84.73367
768.7686
966.3349
−2.90966
−2.90837
−2.90764
−2.86151
−2.79320
−2.73491
−2.33511
−0.58565
1.52752
3.03926
11.20842
37.91648
65.75751
84.60239
768.5391
966.0951
−2.90483
−2.90358
−2.90289
−2.85816
−2.79116
−2.73364
−2.33581
−0.58245
1.53956
3.05817
11.26517
38.08790
66.03580
84.94788
770.4910
968.3932
−2.08231
−2.18772
−2.21064
−2.56811
−2.78140
−2.90948
−3.44367
−4.54168
−5.31616
−5.74199
−7.33773
−10.20901
−12.09497
−13.10498
−27.19100
−29.33698
−2.13734
−2.24267
−2.26552
−2.62156
−2.83493
−2.96349
−3.50238
−4.61775
−5.40773
−5.84283
−7.47698
−10.42916
−12.37538
−13.41981
−28.12546
−30.38698
−2.13855
−2.24357
−2.26631
−2.62180
−2.83575
−2.96467
−3.50386
−4.61440
−5.39922
−5.83143
−7.45563
−10.39149
−12.32654
−13.36467
−27.95970
−30.20117
−2.13508
−2.24045
−2.26330
−2.62091
−2.83658
−2.96675
−3.51208
−4.63229
−5.41799
−5.84823
−7.44960
−10.28980
−12.12807
−13.10402
−26.15762
−28.06336
−2.00274
−2.11480
−2.13603
−2.45423
−2.64415
−2.75888
−4.24997
−4.96603
−5.36068
−9.51569
−11.27323
−12.21479
−25.35462
−27.35677
−2.05764
−2.17116
−2.19274
−2.51315
−2.70305
−2.81783
−4.32273
−5.05093
−5.45305
−9.70756
−11.51580
−12.48636
−26.15001
−28.25025
−2.06192
−2.17478
−2.19618
−2.51430
−2.70317
−2.81753
−4.31868
−5.04363
−5.44360
−9.66842
−11.46230
−12.42486
−25.95352
−28.02904
−2.05915
−2.17102
−2.19218
−2.50754
−2.69480
−2.80795
−4.28822
−5.00052
−5.39247
−9.47577
−11.17646
−12.08100
−24.25357
−26.03943
WFMb
1s 2
1s2p−1
1s3d−2
a
b
−2.861679996
−2.860417861
−2.859709376
−2.814450946
−2.746839677
−2.688884848
−2.289144423
−0.532445132
1.591274097
3.110633781
11.319608967
38.143903320
66.092085756
85.004177725
−2.2353
−2.8301
−5.4000
−12.1011
−26.1264
−2.1659
−2.6871
−4.9866
−11.1540
−24.2283
−2.903351
−2.902083
−2.855859
−2.787556
−2.729508
−2.329780
−0.574877
3.064582
11.267051
38.076320
84.918313
−2.133149
−2.238504
−2.620021
−2.835619
−2.965504
−3.508911
−4.625491
−5.839475
−7.440556
−10.28410
−13.10478
−2.055635
−2.166519
−2.502362
−2.689916
−2.803296
−4.284050
−5.387931
−9.476057
−12.088566
For the 1s 2 state, data are from Ref. [17]; for other states, data are from Ref. [20].
Data are from Ref. [21].
the local and semilocal exchange functionals overestimate
|Ex | when B is very large, thereby overcompensating
for the self-interaction error. (Even for B = 0, spurious
compensation behavior from simple xc functionals has been
known for a long time. LDA xc gets the H atom total energy
roughly right by having a completely spurious nonzero
Ec [63].)
Figure 5 displays the ratio |Ex /ESI | for the fully spinpolarized state of the He atom, 1s2p−1 . Since the antisymmetric wave function makes the two electrons well separated
022504-7
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
TABLE IV. Atomic energies of the lithium atom in B fields by different methods. (Energies in hartree, B in a.u. The chosen xc functionals
for LDA, GGA, and MGGA calculations are PW92, PBE, and TPSS functionals, respectively. For HGGA, we use exact exchange and TPSS
correlation functional.)
HF, lit.a
LDA
GGA
MGGA
HGGA
−7.43275
−7.46856
−7.47740
−7.40878
−7.19621
−6.08810
−5.90114
−3.35784
3.49120
71.80766
939.55235
−7.43275
−7.46857
−7.47741
−7.40879
−7.19621
−6.08811
−5.90113
−3.35777
3.49120
71.807
939.54
−7.34328
−7.37939
−7.39479
−7.33924
−7.14208
−6.04813
−5.86217
−3.32762
3.50491
71.67426
938.1123
−7.46217
−7.49852
−7.51699
−7.46832
−7.28261
−6.21946
−6.03738
−3.54382
3.21595
71.03954
936.0146
−7.48906
−7.52548
−7.54402
−7.49671
−7.31483
−6.26245
−6.08170
−3.60153
3.13744
70.88976
935.8673
−7.48223
−7.51840
−7.52833
−7.45885
−7.24582
−6.13916
−5.95240
−3.41095
3.43600
71.74951
939.4964
−7.477766
−7.517154
−7.528055
−7.458550
−7.244919
−6.136918
−5.949297
−3.406556
−7.4763360 −7.4777957
−7.5122102 −7.5137817
−7.5216127 −7.5235946
−7.4529046
−7.2397460
−7.36507
−7.44174
−7.58787
−7.66652
−7.66245
−6.94229
−6.79515
−4.61775
1.70566
68.17349
930.84309
−7.36509
−7.44176
−7.58790
−7.66653
−7.66246
−6.94230
−6.79517
−4.61777
1.70565
68.1735
930.84308
−7.27909
−7.35647
−7.50718
−7.58789
−7.58026
−6.83784
−6.68788
−4.48455
1.86892
68.30834
929.6385
−7.39752
−7.47501
−7.62468
−7.70476
−7.69981
−6.97522
−6.82790
−4.65417
1.64374
67.80872
927.9366
−7.42270
−7.49961
−7.64670
−7.72652
−7.72379
−7.00375
−6.85670
−4.68388
1.61495
67.78076
927.9620
−7.41506
−7.49225
−7.64052
−7.72182
−7.72297
−7.01341
−6.86726
−4.69801
1.61654
68.06848
930.7310
−7.407126
−7.484773
−7.634547
−7.716679
−7.715709
−7.002346
−6.855410
−4.684076
−7.4088037 −7.4097907
−7.4858382 −7.4869343
−7.6341245 −7.6362483
−7.7151944
−7.7137805
State B (a.u.) HF, present
WFM 1b
WFM 2c
WFM 3d
1s 2 2s
0
0.1
0.5
1
2
5
5.4
10
20
100
1000
1s 2 2p−1
0
0.1
0.5
1
2
5
5.4
10
20
100
1000
1s2p−1 3d−2
0
0.1
1
2
5
10
100
1000
−5.9448544 −5.8555577
−3.4020661
3.4446412
71.7573135
−6.8503299 −6.8361629
-4.6811073
1.6398139
−5.08377 −5.08379 −5.00614 −5.09830 −5.09921 −5.09177 −5.142319
−5.32138 −5.32140 −5.24405 −5.33574 −5.33547 −5.33166 −5.341030
−6.57079 −6.57081 −6.48534 −6.58053 −6.57729 −6.58846 −6.582361
−7.52002 −7.52003 −7.42928 −7.52682 −7.52159 −7.54072 −7.530125
−9.57693 −9.57694 −9.47135 −9.58275 −9.56775 −9.60391 −9.591769
−11.93900 −11.93902 −11.82087 −11.95879 −11.92698 −11.97189 −11.957294
−27.01926 −27.0192 −27.07687 −27.47894 −27.36026 −27.06791
−60.05888 −60.0589 −62.02394 −63.34180 −63.0414 −60.1164
a
Data are from Ref. [16].
Data are from Ref. [24].
c
Data are from Ref. [25].
d
Data are from Ref. [23].
b
in space, the true two-electron exchange contributes only a
small fraction of Ex . About 90% or more of Ex is simply to
counter ESI . Here again, both the PBE and TPSS exchange
functionals do reasonably well when B is small, but fail
badly for larger B. HF exchange exhibits nonmonotonic
behavior in Ex /ESI , with near saturation at B = 10 a.u.
All three approximate xc functionals checked here lack this
characteristic. Instead, they grossly overshoot the X energy
in magnitude at large B. Hence, ESI cannot be expected to
be canceled numerically by a local or semilocal X functional
when B is significant (even if that cancellation were to be from
some mixture of X and C as in the B = 0 case). Because of
the field-induced transverse compression of the density, the
self-interaction problem appears in particularly severe form in
nB-DFT calculations.
Since ExXX , Eq. (10) and EcTPSS both are free from electron
self-interaction by design, the HGGA combination of the
two, Eq. (9), also is self-interaction free. Limitations of this
combination were addressed in Sec. II, but they do not pose a
problem to our present study on magnetized atoms and ions.
HGGA results are presented in Tables I, III, and IV, and in
the Supplemental Material [59]. For Z 4, our HGGA results
agree better with correlated wave-function results [19–27] than
those from other DFT approximations we investigated. For
Z 5, we did not find correlated wave-function results with
which to compare. Table I, however, shows that our HGGA
results are the only ones from DFT which come quite close to
the best available QMC calculations for the atomic energies
of the B and C atoms in large B fields. The comparative
computational cost advantage of HGGA is obvious.
Detailed insight into the excellent behavior of this HGGA
functional comes from comparison with the QMC study in
Ref. [28]. Those authors observed that density-functional
approximations cannot guarantee in all cases to “produce an
022504-8
COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . .
PHYSICAL REVIEW A 90, 022504 (2014)
FIG. 1. (Color online) Atomic total energies for the He atom fully
spin-polarized state 1s2p−1 in magnetic fields from Hartree-Fock and
DFT calculations with different functionals.
FIG. 3. (Color online) One-electron and two-electron energies
for the He atom fully spin-polarized state 1s2p−1 in magnetic
fields from Hartree-Fock and different DFT exchange-correlation
functionals.
upper bound on the ground-state energy in magnetic fields
. . . .” Such failures are examples of the well-known issue
of possible loss of N representability in an approximate xc
functional. However, comparison of the HGGA results with the
released phase diffusion QMC (RPDQMC) results of Ref. [28]
show that Etot,HGGA (B) never is below Etot,RPDQMC (B) (see
Table I). The salient point is that while the key features of
HGGA do not guarantee N representability, they certainly
make it plausible. By construction Ex in HGGA is variationally
exact. While EcTPSS is not guaranteed to be variationally exact,
it is self-interaction-free, hence does not rely upon spurious
overly strong (too negative) C to address SIE. Furthermore,
Ec is a small fraction of ExXX , so one would not expect EcTPSS
to drive the total energy below the exact energy. The numerical
evidence is wholly consistent with this, so operationally there
is evidence that HGGA is N representable.
All of the foregoing DFT calculations rest on the nB-DFT
approximation, with current effects in Exc neglected. The
excellent agreement of our HGGA results with those from the
most sophisticated correlated WFMs and QMC calculations
reveals two messages. One is that the current dependence of
Exc apparently can be neglected safely over a wide range of
experimentally accessible B fields, at least for atomic systems.
The nB-DFT approximation is a useful approach in which the
major B-field effects are included. Another message is that
the development of ground-state DFT functionals with B = 0,
especially the nonlocal HGGA functionals which treat the
self-interaction explicitly has in fact made accurate B > 0 DFT
calculations available, at least for light atoms. This suggests
that nonlocal DFT functionals may be a better starting point
for devising a full-fledged CDFT functional which would also
encompass the current effect in Exc .
FIG. 2. (Color online) Atomic total energies for the Li atom fully
spin-polarized state 1s2p−1 3d−2 in magnetic fields from HartreeFock and DFT calculations with different functionals.
FIG. 4. (Color online) The ratio of exchange energy Ex over
electron self-interaction energy ESI for the He atom singlet state 1s 2
in magnetic fields from Hartree-Fock and different DFT exchange
functionals. The ratio is exactly −1 for Hartree-Fock calculation.
022504-9
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
TABLE V. CDFT corrections to the LDA PW2 results within
VRG approximation (parameters ncutoff = 0.001a0−3 , αcutoff = 2a0−1
VRG
are used for the cutoff function, Exc
in hartree).
Atom
FIG. 5. (Color online) The ratio of exchange energy Ex over
electron self-interaction energy ESI for the He atom fully spinpolarized state 1s2p−1 in magnetic fields from Hartree-Fock and
different DFT exchange functionals.
C. Results for the CDFT VRG approximation
For CDFT, there does not appear to have been any fully
self-consistent calculations on atoms in large B fields based
on the VRG functional [Eq. (11)]. There is related work on
the perturbative implementation of the VRG functional by
Orestes, Marcasso, and Capelle [33] for atomic ionization
energies in vanishing B field. Gross and collaborators compared spin-DFT and CSDFT in several different systems,
including light atoms (for the spurious energy splittings
between current-carrying and zero-current ground states) [47],
magnetic and nonmagnetic solids for the spin-orbit splittings,
and orbital magnetic moments [64], both at B = 0 and the
X-only functional ExXX implemented via the OEP procedure.
They found little differences in the outcome from including
jp dependence or not. For a two-dimensional quantum dot
(2D-QD) in B, they added self-interaction-corrected LSDA
correlation energies post hoc to the bare ExXX results to
compare with QMC energies. Again they found little effect
from the current density [65]. Similar results with the VRG
functional recently were presented by Tellgren et al. [34]
for molecular magnetizabilities, hypermagnetizabilities, and
NMR shielding constants for B up to 1 a.u.
An important formal point is that the VRG functional
VRG
was constructed as an additive term Exc
to provide the
current contribution to Exc at the LDA level of refinement.
That is, VRG is supposed to be used in conjunction with
an LDA functional. The VRG functional also is negative
definite by design. Since LDA energies can be either above
or below accurate results in both low- and high-field regimes
(see tables), the VRG correction cannot always be a proper
correction to bring LDA values closer to exact total energies.
Table V lists representative results for the PW92 + VRG
approximation for the fully spin-polarized states of the He
and Li atoms. An estimate of the effects of jp contributions
is to evaluate the VRG functional using the LDA Kohn-Sham
orbitals. Results from that post-SCF procedure are listed in
VRG
. They are very close to SCF
Table V as non-SCF Exc
State
He
1s2p−1
Li
1s2p−1 3d−2
B (a.u.)
VRG
Non-SCF Exc
VRG
SCF Exc
0
0.5
1
5
10
100
−0.0022
−0.0045
−0.0077
−0.036
−0.074
−0.81
−0.0021
−0.0047
−0.0081
0
2
5
10
−0.0070
−0.027
−0.065
−0.129
−0.0071
−0.029
VRG
Exc
whenever self-consistency can be reached. Evidently
full self-consistency in CDFT calculations has little effect on
the current correction term and, hence, the total energy. Since
SCF convergence problems occur when B is sufficiently large
because of the divergent vorticity ν, we have used non-SCF
VRG
Exc
values in those cases for which self-consistency could
not be reached. Similar observations were reported by Helbig
et al. in 2D-QD, except that they used LSDA for the correlation
energy [65]. This behavior is consistent with the variational
nature of DFT. A small change in the density near the
self-consistent density only results in a second-order energy
change. Presumably the change in density from addition of
current effects is small.
Table V shows that VRG corrections depend strongly
upon particular atomic configurations, but the overall values
are small in magnitude compared to the variation among
HF, LDA, and gradient-dependent DFT results when the B
field is not too large (below roughly 1 a.u.). Within each
configuration, the correction increases with increasing B field.
For B < 5 a.u., the VRG correction is rather small compared
with the discrepancy between the LDA and WFM results and
does not bring them significantly closer. At B = 10 a.u., the
VRG functional fortuitously gives a good correction to the Li
atom 1s2p−1 3d−2 state, but at B = 100 a.u., LDA itself gives
an accurate energy for the He atom 1s2p−1 state. At that point,
the VRG correction can only degrade the results. Overall, there
is no obviously systematic behavior.
As a by-product of numerical issues, we were forced to
study the vorticity ν, which turns out to be a rather difficult
variable to handle. For example, in Fig. 6 we display various
quantities in CDFT along z and radially along ρ for the He
atom 1s2p−1 state in B = 1 a.u. Because jp (0,z) = 0, it is not
displayed. However, ν is not zero along z. Instead, it diverges
at the two poles of the atom and so does the VRG xc energy
density gLCH (n) |ν|2 . Notice that n(r) decays rapidly along z.
At z = 3 bohrs, it is smaller than 10−4 bohr−3 , but the VRG
xc energy density gLCH (n) |ν|2 becomes increasingly large. It
seems rather peculiar that the largest corrections from the VRG
xc functional is at those places where both the electron density
and the current density are nearly zero. Further analysis shows
that this strange behavior is rooted in the choice of ν as a basic
variable in the VRG functional.
022504-10
COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . .
PHYSICAL REVIEW A 90, 022504 (2014)
we view the CDFT corrections for total energies as subject to
caution; they are semiquantitative at best.
V. DISCUSSION AND CONCLUSIONS
FIG. 6. (Color online) Various quantities (electron density n,
paramagnetic current density jp , vorticity ν, and the magnitude of the
current correction to the exchange-correlation energy density, |gν 2 |,
in the VRG functional) for the He atom 1s2p−1 state in B = 1 a.u.
All quantities are evaluated from the LDA KS orbitals and plotted
along the z and ρ axes (cylindrical coordinates).
One may argue that at the places of such a low density the
prefactor g(n) is not well defined in practice because gLCH (n)
(or other fitted functions) only used data points with rs 10
bohrs [49]. This is indeed a problem, but unfortunately we do
not know the form of g(n) in the low-density limit. The picture
of HEG at a very low electron density in a strong B field may
be inappropriate anyway.
To address this divergence problem numerically, we introduced a rapidly decaying cutoff function gcutoff . Details
were given in the methodology section. However, there is
no obvious, physically grounded set of criteria for choosing
the cutoff parameters ncutoff and αcutoff . We tested different
parameter sets by using the He atom 1s2p−1 state in B = 1 a.u.
Different cutoff parameters generate quite different results
see Table VI). This outcome is distinct from the findings
in Ref. [34] presumably for two reasons. Our study focuses
on total energies, theirs on response (e.g., magnetizability).
Second, we have gone to much higher fields. The test results
shown in Table VI are at B = 1 a.u., the highest field used in
Ref. [34] but at the low end of the range we have examined.
Comparison with the He 1s2p−1 LDA energy at B = 1 a.u.
(−2.909 48 hartree; recall Table III) shows the undesirable
effect upon total energies of the dependence of the VRG
correction upon regularization. Because of that dependence,
TABLE VI. Effect of cutoff parameters on CDFT VRG corrections for the He atom 1s2p−1 state in magnetic field B = 1 a.u.
(energies in hartree).
ncutoff (units of a0−3 )
0.005
0.001
0.001
0.0001
0.00001
αcutoff (units of a0−1 )
VRG
Non-SCF Exc
2.0
2.0
1.0
2.0
2.0
−0.004
−0.008
−0.010
−0.025
−0.064
It is perhaps unsurprising that the VRG functional fails
when applied to atomic systems in large B. The functional
was developed from the study of the moderately dense to
dense HEG in low B, for which Landau orbitals were used
as approximations. This physical situation is quite different
from a finite system such as an atom. First, n(r) and jp (r)
vary considerably within an atom, and the low-density regions
(rs > 10 a0 ) are non-negligible. Secondly, there is not a direct
relationship between jp (r) and the external B field as there
is for the HEG. The question whether the electron gas
remains homogeneous after imposing a substantial B field
is even unclear. If the field were to induce some form of
crystallization, the basic picture on which the VRG functional
is based is lost. As Pittalis et al. pointed out, the B field
can induce derivative discontinuities in the xc energy density.
Such discontinuities make practical calculations extremely
difficult [47]. Our analysis and numerical studies also suggest
the picture of Landau orbitals used for the HEG may not
be applicable at all for the atomiclike systems. Unlike the
LDA, also based on the HEG, the VRG functional may be
too simple to encompass the essential physics of the atomic
systems. Alternatives are admittedly difficult to imagine.
A more fundamental question is whether ν is a suitable
basic variable in gauge-invariant CDFT as Vignale and Rasolt
required [6]. While appealing on purely theoretical grounds,
our numerical results on atomic systems in B suggest ν is,
at best, an awkward choice. It is ν as the basic variable in
the VRG functional which drives unphysical results and the
need for cutoffs. While ν appears to be a misbehaved quantity
for atomic systems, jp is still a well-behaved variable, but
it suffers from gauge dependence. Recently, Tellgren et al.
reexamined the possibility of using (n,jp ) as basic variables,
and found that the gauge problem can be circumvented by
using modified conjugated potentials [66]. Perhaps one should
take a step back from ν to jp as a basic variable to conceive
viable CDFT functionals. We note, however, that this is not the
route suggested in Ref. [67]. Pan and Sahni also formulated
CDFT by using the electron density n and physical current
density j [68], providing another interesting idea, though with
some restrictions [69], the consequences of which remain to
be explored (see, for example, Ref. [70]).
There is a larger issue. To a considerable extent, the present
study calls into question the utility of CDFT (in any form)
itself. Exact exchange combined with the TPSS correlation
functional yields a self-interaction-free HGGA functional,
Eq. (9). Energies from this functional used as an nB-DFT
approximation agree quite well with results from correlated
WFMs and QMC. That agreement indicates that the current
correction to Exc from CDFT is rather small over a wide range
of B. As a result, the nonlocal HGGA functional may be both
a practical alternative to CDFT (at least over the range of fields
which we have considered) as well as a better starting point
than LDA to construct full-fledged CDFT functionals.
One final remark is pertinent to the results presented here.
Neither relativistic effects nor finite nuclear mass effects
022504-11
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
have been considered. Those effects can be important for
matter in superstrong fields (B > 104 a.u.), in which regime
the adiabatic approximation will be applicable. But for the
field strengths covered in this paper, both effects should be
negligible.
ACKNOWLEDGMENTS
Helpful discussions with Gao Xianlong (Zhejiang Normal
University), Wenjian Liu (Peking University), X.-Y. Pan
(Ningbo University), and Jian Wang (Huzhou University) are
acknowledged with thanks. This work was funded by Zhejiang
Provincial Natural Science Foundation of China under Grant
No. LY13A050002 (W.Z. and L.Z.), in part by the National
Natural Science Foundation of China Grant No. 11274085
(W.Z.), and in part by U.S. Department of Energy Grant No.
DE-SC-0002139 (S.B.T.). In early portions of this work both
W.Z. and S.B.T. were supported by U.S. National Science
Foundation Grant No. DMR-0218957.
APPENDIX
The KL basis sets start with a base sequence of orbital
exponents defined as those nearest to spherical, i.e., closest to
αj = βj ∀j . The derived sequences compress the transverse
(radial) functions (the αj ) relative to the corresponding
longitudinal (axial) ones (βj ). Thus, the KL basis sets have
several (one to five) exponent sequences of the form
μ2,3 = ±0.2,
μ4,5 = ±0.4,
(A1)
with αj,1 ≈ βj ∀j and KL a parameter to be optimized at one
fixed B. Note that in each sequence the differences between the
transverse (αj ) and longitudinal exponents (βj ) are the same
for all members of that sequence; neither μ nor KL is indexed
by j . To keep the basis size within reason, the second and third
sequences ( = 2,3) are half the length (with doubled spacing)
of the first sequence, while the fourth and fifth sequences
( = 4,5) are one-quarter as long (with quadrupled spacing)
as the base sequence.
We reworked the KL scheme as follows. First was some
numerical experimentation on the H atom. At B = 0 a.u., an
even-tempered Gaussian sequence of length Nb = 16 and rule
of formation [71],
αj = βj ≡ pq j , j = 1,2, . . . ,Nb ,
ln p = a ln(q − 1) + a ,
ln(ln q) = b ln Nb + b ,
a = 0.3243,
b = −0.4250,
:= βj + (1 + μ1 )(B),
(A2)
a = −3.6920,
b = 0.9280,
gives Etot = −0.499 999 92 hartree. However, at B = 10 a.u.,
it gives a 24% error compared to the known value Etot (B =
10) = −1.747 797 163 714 hartree [60]. The task therefore
is to discover as near-optimum asymmetry between the
transverse and longitudinal manifolds as feasible.
Within only the basic KL sequence and the Nb = 16
spherical set, we searched in the parameter space αj to
minimize the total energy of the hydrogen atom at B = 10 a.u.
(A3)
where
b(γ ) = −0.16[tan−1 (γ )]2 + 0.77 tan−1 (γ ) + 0.74.
(A4)
is the reduced field strength for an
In general γ =
effective nuclear charge Zeff . For the innermost electrons, Zeff
is close to the bare nuclear charge. For other electrons, it is
close to the nuclear charge minus the number of inner shell
electrons (screening effects). In practice, however, we found
that nominal Zeff values (bare nuclear charge), such as just
given, are good enough. The expressions for αj , Eq. (A3), and
b(γ ), Eq. (A4), came from fits to the optimized energies of the
H isoelectronic series (H, He+ , Li++ , Be+++ , C5+ , and O7+ in
reduced fields γ = 0.1, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500,
800, 1000, 2000, and 4000. Those optimizations also started
from the Nb = 16 basis set.
Equation (A3) with μ1 = 0 defines the first sequence.
Other sequences are defined analogously with the KL basis
set, μ = ±0.2, ± 0.4 for the second, third, fourth, and fifth
sequences (with increased spacing as mentioned in the main
text), respectively.
Similar to the KL basis sets, we used ETG sequences for the
longitudinal exponents βj ; recall Eq. (A2). Since the magnetic
field does not change the confinement along the z direction, we
expect it to have little effect upon ETG exponents. In cases for
which the electron density has diffuse, nonzero orbital angular
momentum contributions, extrapolation of the tempering to
include a small number of negative j values is used.
Testing confirmed that the resulting basis sets work as
well for multielectron atoms as for one-electron systems. Our
criterion for basis set error is that it be less than 1 μHartree
in the whole B range investigated (0 B 2000 a.u.) for the
hydrogen atom, with respect to more accurate algebraic results
[60]. We find that the HF energies for multielectron atoms from
our basis sets are nearly indistinguishable from full numerical
two-dimensional mesh results [16].
2
B/Zeff
αj, = βj + (1 + μ )KL ,
μ1 = 0.0,
and found that the optimized difference αj − βj is not constant
at fixed B as the KL sequences suggest. Particularly for
the smaller αj values, the fractional difference (αj − βj )/αj
is quite large, as much as ≈ 0.95. The behavior is quite
understandable. A small-exponent basis function samples
a large volume far from the nucleus, where the B field
overpowers the nuclear attraction such that the distortion from
the field-free spherical symmetry will be relatively larger than
for the region sampled by larger-exponent basis functions. In
the limit βj → 0, which can be thought of either as the large
B limit or zero nuclear charge, the electron wave function is a
Landau orbital, with an exponential parameter αj ≡ aB = B/4
(aB is the scale parameter in the Landau orbital). The opposite
limit, βj → ∞, corresponds to B = 0 a.u., for which spherical
symmetry is restored, αj = βj . Thus, a convenient choice for
orbital exponent asphericity is a scaling of aB .
By consideration of such large and small B-field limits and
use of a nonlinear fit to calculated one-electron system results,
we reached a prescription for nearly optimized exponents,
namely,
B
4 βj −1/2
4 βj −2
αj = βj +
+ 1+
4 1+
20
b(γ ) B
b(γ ) B
022504-12
COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . .
PHYSICAL REVIEW A 90, 022504 (2014)
[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964); W.
Kohn and L. J. Sham, ibid. 140, A1133 (1965).
[2] P. B. Jones, Phys. Rev. Lett. 55, 1338 (1985).
[3] B. M. Relovsky and H. Ruder, Phys. Rev. A 53, 4068
(1996).
[4] M. Braun, Phys. Rev. A 65, 033415 (2002).
[5] Z. Medin and D. Lai, Phys. Rev. A 74, 062507 (2006).
[6] G. Vignale and M. Rasolt, Phys. Rev. Lett. 59, 2360 (1987);
,Phys. Rev. B 37, 10685 (1988); G. Vignale, M. Rasolt,
and D. J. W. Geldart, Adv. Quantum Chem. 21, 235 (1990).
[7] Wuming Zhu and S. B. Trickey, J. Chem. Phys. 125, 094317
(2006).
[8] M. Taut, P. Machon, and H. Eschrig, Phys. Rev. A 80, 022517
(2009).
[9] E. H. Lieb and R. Schrader, Phys. Rev. A 88, 032516
(2013).
[10] P. Pr¨oschel, W. R¨osner, G. Wunner, H. Ruder, and H. Herold, J.
Phys. B 15, 1959 (1982).
[11] D. Neuhauser, S. E. Koonin, and K. Langanke, Phys. Rev. A 36,
4163 (1987).
[12] D. Lai, E. E. Salpeter, and S. L. Shapiro, Phys. Rev. A 45, 4832
(1992).
[13] H. Ruder, G. Wunner, H. Herold, and F. Geyer, Atoms in Strong
Magnetic Fields (Springer-Verlag, Berlin, 1994).
[14] J. S. Heyl and L. Hernquist, Phys. Rev. A 58, 3567 (1998); A.
Thirumalai and J. S. Heyl, ibid. 79, 012514 (2009).
[15] M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys. Rev. A 54,
219 (1996).
[16] M. V. Ivanov and P. Schmelcher, Phys. Rev. A 57, 3793 (1998);
,60, 3558 (1999); ,61, 022505 (2000); ,Eur. Phys. J. D 14, 279
(2001); ,J. Phys. B 34, 2031 (2001).
[17] L. B. Zhao and P. C. Stancil, Phys. Rev. A 77, 035401 (2008).
[18] D. Engel and G. Wunner, Phys. Rev. A 78, 032515 (2008).
[19] A. Scrinzi, Phys. Rev. A 58, 3879 (1998).
[20] M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys. Rev. A 59,
2875 (1999).
[21] W. Becken, P. Schmelcher, and F. K. Diakonos, J. Phys. B 32,
1557 (1999); W. Becken and P. Schmelcher, ibid. 33, 545 (2000);
,Phys. Rev. A 63, 053412 (2001).
[22] X. Wang, J. Zhao, and H. Qiao, Phys. Rev. A 80, 053425 (2009).
[23] X. Guan and B. Li, Phys. Rev. A 63, 043413 (2001).
[24] O.-A. Al-Hujaj and P. Schmelcher, Phys. Rev. A 70, 033411
(2004).
[25] X. Wang and H. Qiao, Phys. Rev. A 75, 033421 (2007).
[26] X. Guan, B. Li, and K. T. Taylor, J. Phys. B 36, 2465
(2003).
[27] O.-A. Al-Hujaj and P. Schmelcher, Phys. Rev. A 70, 023411
(2004).
[28] S. B¨ucheler, D. Engel, J. Main, and G. Wunner, Phys. Rev. A
76, 032501 (2007).
[29] C. Schimeczek, S. Boblest, D. Meyer, and G. Wunner, Phys.
Rev. A 88, 012509 (2013).
[30] S. M. Colwell and N. C. Handy, Chem. Phys. Lett. 217, 271
(1994); A. M. Lee, S. M. Colwell, and N. C. Handy, ibid. 229,
225 (1994).
[31] A. M. Lee, N. C. Handy, and S. M. Colwell, J. Chem. Phys. 103,
10095 (1995).
[32] S. M. Colwell, N. C. Handy, and A. M. Lee, Phys. Rev. A 53,
1316 (1996); A. G. Ioannou, S. M. Colwell, and R. D. Amos,
Chem. Phys. Lett. 278, 278 (1997).
[33] E. Orestes, T. Marcasso, and K. Capelle, Phys. Rev. A 68,
022105 (2003).
[34] E. I. Tellgren, A. M. Teale, J. W. Furness, K. K. Lange, U.
Ekstr¨om, and T. Helgaker, J. Chem. Phys. 140, 034101 (2014).
[35] B. G. Johnson, P. M. W. Gill, and J. A. Pople, J. Chem. Phys.
98, 5612 (1993).
[36] I. H. Lee and R. M. Martin, Phys. Rev. B 56, 7197 (1997).
[37] S. Kurth, J. P. Perdew, and P. Blaha, Int. J. Quantum Chem. 75,
889 (1999).
[38] J. P. Perdew and K. Schmidt, Density Functional Theory and
Its Application to Materials, edited by V. Van Doren, C. Van
Alsenoy, and P. Geerlings, AIP. Conf. Proc. No. 577 (AIP,
Melville, NY, 2001), p. 1.
[39] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
[40] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77,
3865 (1996); ,78, 1396 (1997).
[41] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys.
Rev. Lett. 91, 146401 (2003).
[42] O. A. Vydrov, G. E. Scuseria, J. P. Perdew, A. Ruzsinszky, and
G. I. Csonka, J. Chem. Phys. 124, 094108 (2006).
[43] A. G¨orling, Phys. Rev. Lett. 83, 5459 (1999).
[44] S. K¨ummel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008).
[45] L. A. Constantin, J. M. Pitarke, J. F. Dobson, A. Garcia-Lekue,
and J. P. Perdew, Phys. Rev. Lett. 100, 036401 (2008).
[46] J. D. Talman and W. F. Shadwick, Phys. Rev. A 14, 36 (1976).
[47] S. Pittalis, S. Kurth, N. Helbig, and E. K. U. Gross, Phys. Rev.
A 74, 062511 (2006).
[48] J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 45, 101
(1992).
[49] G. Vignale, M. Rasolt, and D. J. W. Geldart, Phys. Rev. B 37,
2502 (1988).
[50] J. Tao and J. P. Perdew, Phys. Rev. Lett. 95, 196403 (2005).
[51] J. Tao and G. Vignale, Phys. Rev. B 74, 193108 (2006).
[52] K. Higuchi and M. Higuchi, Phys. Rev. B 74, 195122 (2006).
[53] M. D. Jones, G. Ortiz, and D. M. Ceperley, Int. J. Quantum
Chem. 64, 523 (1997).
[54] D. Porezag and M. R. Pederson, Phys. Rev. A 60, 2840 (1999).
[55] C. Aldrich and R. L. Greene, Phys. Status Solidi B 93, 343
(1979).
[56] P. Schmelcher and L. S. Cederbaum, Phys. Rev. A 37, 672
(1988).
[57] U. Kappes and P. Schmelcher, J. Chem. Phys. 100, 2878 (1994).
[58] Yu. P. Kravchenko and M. A. Liberman, Int. J. Quantum Chem.
64, 513 (1997).
[59] See Supplemental Material at http://link.aps.org/supplemental/
10.1103/PhysRevA.90.022504 for extensive tabulations of HF
and DFT calculated data for the He, Li, Be, B, and C atoms and
the positive ions, Li+ , Be+ , and B+ .
[60] Y. P. Kravchenko, M. A. Liberman, and B. Johansson, Phys.
Rev. A 54, 287 (1996).
[61] R. W. Danz and M. L. Glasser, Phys. Rev. B 4, 94 (1971).
[62] Weitao Yang, P. Mori-Sanchez, and A. J. Cohen, J. Chem. Phys.
139, 104114 (2013).
[63] S. B. Trickey, Phys. Rev. Lett. 56, 881 (1986).
[64] S. Sharma, S. Pittalis, S. Kurth, S. Shallcross, J. K. Dewhurst,
and E. K. U. Gross, Phys. Rev. B 76, 100401(R) (2007).
[65] N. Helbig, S. Kurth, S. Pittalis, E. R¨as¨anen, and E. K. U. Gross,
Phys. Rev. B 77, 245106 (2008).
[66] E. I. Tellgren, S. Kvaal, E. Sagvolden, U. Ekstr¨om, A. M. Teale,
and T. Helgaker, Phys. Rev. A 86, 062506 (2012).
022504-13
WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY
PHYSICAL REVIEW A 90, 022504 (2014)
[67] K. Higuchi, M. Miyasita, and M. Higuchi, Int. J. Quantum Chem.
110, 2286 (2010).
[68] X.-Y. Pan and V. Sahni, Int. J. Quantum Chem. 110, 2833
(2010); ,J. Phys. Chem. Sol. 73, 630 (2012); ,Phys. Rev. A
86, 042502 (2012); V. Sahni and X.-Y. Pan, ibid. 85, 052502
(2012).
[69] X.-Y. Pan and V. Sahni, Int. J. Quantum Chem. 114, 233
(2014).
[70] A. Laestadius and M. Benedicks, Int. J. Quantum Chem. 114,
782 (2014).
[71] M. W. Schmidt and K. Ruedenberg, J. Chem. Phys. 71, 3951
(1979).
022504-14