Exchange parameters from approximate self-interaction correction scheme
Exchange parameters from approximate self-interaction correction scheme
A. Akande and S. Sanvito
School of Physics and CRANN, Trinity College, Dublin 2, Ireland
(Dated: November 28, 2018)
The approximate atomic self-interaction corrections (ASIC) method to density functional theory is
put to the test by calculating the exchange interaction for a number of prototypical materials, critical
to local exchange and correlation functionals. ASIC total energy calculations are mapped onto an
Heisenberg pair-wise interaction and the exchange constants J are compared to those obtained
with other methods. In general the ASIC scheme drastically improves the bandstructure, which
for almost all the cases investigated resemble closely available photo-emission data. In contrast the
results for the exchange parameters are less satisfactory. Although ASIC performs reasonably well
for systems where the magnetism originates from half-filled bands, it suffers from similar problems
than those of LDA for other situations. In particular the exchange constants are still overestimated.
This reflects a subtle interplay between exchange and correlation energy, not captured by the ASIC.
PACS numbers:
I. INTRODUCTION
Theoretical studies based on density functional the-
ory (DFT) [1, 2] have given remarkable insights into
the electronic and magnetic properties of both molecules
and solids [3]. In particular, a number of these studies
attempt to quantitatively describe the magnetic inter-
action in a broad range of systems including transition
metals [4], hypothetical atomic chains [5, 6], ionic solids
[7, 8, 9], transition metal oxides [10, 11] and transition
metal polynuclear complexes [12, 13, 14]. DFT uses an
effective single-particle picture where spin symmetry is
generally broken. For this reason exchange parameters J
are conventionally extracted by using a mapping proce-
dure, where total energy calculations are fitted to a classi-
cal Heisenberg Hamiltonian [4, 15]. This is then used for
evaluating the Curie or Néel temperatures, the magnetic
susceptibility and for interpreting neutron diffraction ex-
periments [16].
Notably, the accuracy and reliability of the numerical
values of the J ’s depend on the functional used for the ex-
change and correlation (XC) energy, being this the only
approximated part of the DFT total energy [17]. Cal-
culations based on well-known local functionals, namely
the local density approximation (LDA) and the gener-
alised gradient approximation (GGA), are successful with
itinerant magnetism in transition metals [4], but largely
over-estimates the Heisenberg exchange parameters in
many other situations [7, 8, 9, 11, 14]. Additional cor-
rections based on the kinetic energy density (metaGGA)
[18] marginally improves the agreement with experiments
[9], although an extensive investigation over several solid
state systems has not been carried out so far.
These failures are usually connected to the local char-
acter of the LDA, which is only weakly modified by con-
structing XC potentials including the gradient, or higher
derivative of the charge density. A direct consequence is
that the charge density is artificially delocalized in space,
leading to an erroneous alignment of the magnetic bands.
These are also artificially broadened. A typical example
is that of NiO, which LDA predicts as Mott-Hubbard in-
stead of charge-transfer insulator. Clearly a qualitative
failure in describing the ground state results in an erro-
neous prediction of the exchange parameters.
One of the reasons behind the inability of LDA and
GGA of describing localized charge densities is attributed
to the presence of the self-interaction error (SIE) [19].
This originates from the spurious Coulomb interaction of
an electron with itself, which is inherent to local func-
tionals. Hartree-Fock (HF) methods, in the unrestricted
or spin polarised form, are SIE free and produce sys-
tematic improved J parameters. However, these meth-
ods lack of correlation and usually overcorrect. A typ-
ical example is the antiferromagnetic insulator KNiF3
for which HF predicts a nearest neighbour J of around
2 meV [7, 20, 21, 22, 23] against an experimental value
of 8.6 meV [24]. Direct SIE subtraction, convention-
ally called self-interaction corrected (SIC) LDA, also im-
proves the results and seems to be less affected by over-
correction [5, 25]. Similarly hybrid-functionals, which
mix portions of HF exchange with the local density ap-
proximation of DFT, perform better than local function-
als and in several situations return values for J in close
agreement with experiments [7, 8].
It is important to note that both methods based non-
local exchange or SIC, are computationally demanding
and thus their application to the solid state remains
rather limited. It is then crucial to develop practical com-
putational schemes able to provide a good estimate of the
exchange parameters for those systems critical to LDA,
which at the same time are not numerically intensive.
Based on the idea that most of the SIE originates from
highly localized states, with a charge distribution resem-
bling those of free atoms, Vogel et al. [26] proposed a sim-
ple SIC scheme where the corrections are approximated
by a simple on-site term. This method was then gener-
alized to fractional occupation by Filippetti and Spaldin
[27] and then implemented in a localized atomic orbital
code for large scaling by Pemmaraju et al. [28]. Despite
its simplicity the method has been successfully applied to
a number of interesting physical systems including , tran-
sition metal monoxides [27, 29], silver halides [30], no-
ble metal oxides [31], ferroelectric materials [27, 32, 33],
high-k materials [34], diluted magnetic semiconductors
[35, 36] and also to quantum transport [37, 38].
The method is strictly speaking not variational, in the
sense that a functional generating the ASIC potential via
variational principle is not available. However, since typ-
ically the LDA energy is a good approximation of the
exact DFT energy, although the LDA potential is rather
different from the exact KS potential, a “practical” def-
inition of total energy can be provided. In this work we
evaluate the ability of this approximated energy in de-
scribing exchange parameters for a variety of magnetic
systems.
II. THE ATOMIC SIC METHOD
The seminal work of Perdew and Zunger [19] pioneered
the modern theory of SIC. The main idea is that of sub-
tracting directly the spurious SI for each Kohn-Sham
(KS) orbital ψn. The SIC-LDA [39] XC energy thus
writes
ESICxc [ρ
↑, ρ↓] = ELDAxc [ρ
↑, ρ↓]−
occupied∑
δSICn , (1)
where ELDAxc [ρ
↑, ρ↓] is the LDA-XC energy and δSICn is the
sum of the self-Hartree and self-XC energy associated to
the charge density ρσn = |ψσn|2 of the fully occupied KS
orbital ψσn
δSICn = U [ρ
n] + E
xc [ρ
n, 0] . (2)
Here U is the Hartree energy and σ is the spin index.
The search for the energy minimum is not trivial, since
ESICxc is not invariant under unitary rotations of the occu-
pied KS orbitals. As a consequence the KS method be-
comes either non-orthogonal or size-inconsistent. These
problems however can be avoided [40, 41, 42] by intro-
ducing a second set of orbitals φσn related to the canonical
KS orbitals by a unitary transformation M
ψσn =
Mσnmφ
m . (3)
The functional can then be minimized by varying both
the orbitals ψ and the unitary transformationM, leading
to a system of equations
n = (H
0 + ∆v
n = �
σ,SIC
n , (4)
ψσn =
Mnmφσm , (5)
∆vSICn =
MnmvSICm
vSICm P̂
m , (6)
where Hσ0 is the LDA Hamiltonian, P̂
n(r) =
φσm(r)〈φσm|ψσn〉 and vSICn = u([ρn]; r) + vσ,LDAxc ([ρ↑n, 0]; r),
with u and vσ,LDAxc the Hatree and LDA-XC potential
respectively.
In equation (4) we have used the fact that at the en-
ergy minimum the matrix of SIC KS-eigenvalues �σ,SICnm
is diagonalized by the KS orbitals ψn. Importantly such
minimization scheme can be readily applied to extended
systems, without loosing the Bloch representation of the
KS orbitals [43, 44].
The ASIC method consists in taking two drastic ap-
proximations in equation (4). First we assume that the
orbitals φm, that minimize the SIC functional are atomic-
like orbitals Φσm (ASIC orbitals) thus∑
vSICm (r)P̂
m → α
ṽσSICm (r)P̂
m , (7)
where ṽσSICm (r) and P̂
m are the SIC potential and the
projector associated to the atomic orbital Φσm. Secondly
we replace the non-local projector P̂Φm with its expecta-
tion value in such a way that the final ASIC potential
reads
vσASIC(r) = α
ṽσSICm (r)p
m , (8)
where pσm is the orbital occupation (essentially the spin-
resolved Mülliken orbital population) of Φm.
Note that in the final expression for the potential a
factor α appears. This is an empirical scaling term that
accounts for the fact that the ASIC orbital Φ in general
do not coincide with those that minimize the SIC func-
tional (1). By construction α = 1 in the single particle
limit, while it vanishes for the homogeneous electron gas.
Although in general 0 < α < 1, extensive testing [28]
demonstrates that a value around 1 describes well ionic
solids and molecules, while a value around 1/2 is enough
for mid- to wide-gap insulators. In the following we will
label with ASIC1/2 and ASIC1 calculations obtained re-
spectively with α = 1/2 and α = 1.
Finally we make a few comments over the total energy.
As pointed out in the introduction the present theory is
not variational since the KS potential cannot be related
to a functional by a variational principle. However, since
typical LDA energies are more accurate than their corre-
sponding KS potentials, we use the expression of equation
(1) as suitable energy. In this case the orbital densities
entering the SIC are those given by the ASIC orbital Φ.
Moreover, in presenting the data, we will distinguish re-
sults obtained by using the SIC energy (1) from those
obtained simply from the LDA functional evaluated at
the ASIC density, i.e. without including the δn correc-
tions (2).
III. RESULTS
All our results have been obtained with an implemen-
tation of the ASIC method [28] based on the DFT code
Siesta [45]. Siesta is an advanced DFT code using pseu-
dopotentials and an efficient numerical atomic orbital ba-
sis set. In order to compare the exchange parameters ob-
tained with different XC functionals we consider the LDA
parameterization of Ceperly and Alder [46], the GGA
functional obtained by combining Becke exchange [47]
with Lee-Yang-Parr correlation [48] (BLYP), the nonem-
pirical Purdew, Burke and Ernzerhof (PBE) GGA [49],
and the ASIC scheme as implemented in reference [28].
Calculations are performed for different systems crit-
ical to LDA and GGA, ranging from molecules to ex-
tended solids. These include hypothetical H-He atomic
chains, the ionic solid KNiF3 and the transition metal
monoxides MnO and NiO. DFT total energy calculations
are mapped onto an effective pairwise Heisenberg Hamil-
tonian
HH = −
Jnm~Sn · ~Sm , (9)
where the sums runs over all the possible pairs of spins.
In doing this we wish to stress that the mapping is a con-
venient way of comparing total energies of different mag-
netic configurations calculated with different function-
als. In this spirit the controversy around using the spin-
projected (Heisenberg mapping) or the non-projected
scheme is immaterial [5, 50, 51].
A. H-He chain
As an example of molecular systems, we consider H-
He monoatomic chains at a inter-atomic separation of
1.625 Å (see figure 1). This is an important benchmark
for DFT since the wave-function is expected to be rather
localized and therefore to be badly described by local
XC functionals. In addition the system is simple enough
to be accessible by accurate quantum chemistry calcula-
tions.
As basis set we use two radial functions (double-ζ)
for the s and p angular momenta of both H and He,
while the density of the real-space grid converges the
self-consistent calculation at 300 Ry. Here we consider
all possible Heisenberg parameters. Thus the triangular
molecule (Fig.1a) has only one nearest neighbour param-
eter Ja12, the 5-atom chain (Fig.1b) has both first J
12 and
second neighbour Jb13 parameters, and the 7-atom chain
(Fig.1c) has three parameters describing respectively the
nearest neighbour interaction with peripheral atoms Jc12,
the nearest neighbour interaction between the two middle
atoms Jc23 and the second neighbour interaction J
Following reference [5], accurate calculations based
on second-order perturbation theory (CASPT2) [12] are
used as comparison. The quality of each particular func-
tionals is measured as the relative mean deviation of the
nearest neighbour exchange parameters only (Ja12, J
FIG. 1: (Color on-line) H-He-H chains at an inter-atomic dis-
tance of 1.625Å.
Method Ja12 J
13 δ (%)
CASPT2 -24 -74 -0.7 -74 -79 -0.7 0
SIC-B3LYP -31 -83 -0.2 -83 -88 -0.3 16
LDA -68 -232 -6 -234 -260 -6 210
PBE -60 -190 -1.8 -190 -194 -1.6 152
BLYP -62 -186 -2 -186 -200 -1 147
ASIC1 -36 -112 -1 -110 -122 -0.6 51
ASIC1/2 -44 -152 -1 -152 -168 -1.4 101
ASIC∗1 -40 -128 -0.6 -128 -142 -1.0 73
ASIC∗1/2 -50 -170 -1.4 -170 -190 -1.8 127
TABLE I: Calculated J values (in meV) for the three different
H–He chains shown in Fig.1. The CASPT2 values are from
reference [12], while the SIC-B3LYP are from reference [5].
The last two rows correspond to J values obtained from the
LDA energy calculated at the ASIC density.
Jc12, J
23), since those are the largest ones
|Ji − JCASPT2i |
|JCASPT2i |
. (10)
Our calculated J values and their relative δ are pre-
sented in table I, where we also include results for a fully
self-consistent SIC calculation over the B3LYP functional
(SIC-B3LYP) [5]. It comes without big surprise that the
LDA systematically overestimates all the exchange pa-
rameters with errors up to a factor 6 for the smaller J
(Jb13 and J
13) and an average error δ for the largest J
of about 200%. Standard GGA corrections considerably
improve the description although the J ’s are still system-
atically larger than those obtained with CASPT2. Note
that the results seem rather independent of the particular
GGA parameterization, with PBE and BLYB producing
similar exchange constants. This is in good agreement
with previous calculations [5].
SIC in general dramatically improves the LDA and
GGA description and our results for ASIC1 are reason-
ably close to those obtained with the full self-consistent
procedure (SIC-B3LYP). This is an interesting result,
considering that our ASIC starts from a local exchange
functional, while B3LYP already contains non-local con-
tributions. We also evaluate the J parameters by using
the LDA energy evaluated at the ASIC density (last two
rows in table I). In general this procedure gives J ’s larger
than those obtained by using the energy of equation (1),
meaning that the δSICn contributions reduce the J values.
It is then clear that the ASIC scheme systematically
improves the J values as compared to local functionals.
The agreement however is not as good as the one ob-
tained by using a fully self-consistent SIC scheme, mean-
ing that for this molecular system the ASIC orbitals are
probably still not localized enough. This can alter the
actual contribution of δSICn to the total energy and there-
fore the exchange parameters.
B. Ionic antiferromagnets: KNiF3
Motivated by the substantial improvement of ASIC
over LDA, we then investigate its performances for real
solid-state systems, starting from KNiF3. This is a pro-
totypical Heisenberg antiferromagnet with strong ionic
character, a material for which our ASIC approxima-
tion is expected to work rather well [28]. It is also a
well studied material, both experimentally [24, 52] and
theoretically [7, 9, 21, 22], allowing us extensive com-
parisons. The KNiF3 has cubic perovskite-like structure
with the nickel atoms at the edges of the cube, flourine
atoms at the sides and potassium atoms at the center (see
Fig.2). At low temperature, KNiF3 is a type II antiferro-
magnetic insulator consisting of ferromagnetic (111) Ni
planes aligned antiparallel to each other. For our calcu-
lations we use a double-ζ polarized basis for the s and p
orbitals of K, Ni and F, a double-ζ for the 3d of K and Ni,
and a single-ζ for the 3d of F. Finally, we use 5×5×5 k-
points in the full Brillouin zone and the real-space mesh
cutoff is 550 Ry. Note that the configuration used to
generate the pseudopotential is that of Ni2+, 4s13d7.
We first consider the band-structure as obtained with
LDA and ASIC. For comparison we also include results
obtained with LDA+U [53, 54] as implemented in Siesta
[55]. In this case we correct only the Ni d shell and we
fix the Hubbard-U and Hund’s exchange-J parameters by
fitting the experimental lattice constant (a0 = 4.014 Å).
The calculated values are U=8 eV and J=1 eV. The
bands obtained with the three methods and the corre-
sponding orbital projected density of states (DOS) are
presented in figures 3 and 4 respectively.
All the three functionals describe KNiF3 as an insula-
tor with bandgaps respectively of 1.68 eV (LDA), 3.19 eV
(ASIC1), and 5.0 eV (LDA+U). An experimental value
for the gap is not available and therefore a comparison
cannot be made. In the case of LDA and ASIC the gap
is formed between Ni states, with conductance band bot-
FIG. 2: (Color on-line) Cubic perovskite structure of KNiF3.
Color code: blue=Ni, red=F, Green=K.
FIG. 3: Band structure for type II antiferromagnetic KNiF3
obtained with a) LDA, b) ASIC1 and c) LDA+U (U=8 eV
and J=1 eV). The valence band top is aligned at E=EF=0 eV
(horizontal line).
tom well described by eg orbitals. These are progressively
moved upwards in energy by the SIC, but still occupy the
gap. Such feature is modified by LDA+U which pushes
the unoccupied eg states above the conductance band
minimum, which is now dominated by K 4s orbitals.
In more detail the valence band is characterized by a
low-lying K 3p band and by a mixed Ni-3d/F 2p. While
the K 3p band is extremely localized and does not present
substantial additional orbital components the amount of
mixing and the broadening of the Ni-3d/F 2p varies with
the functionals used. In particular both LDA and ASIC
predict that the Ni 3d component occupies the high en-
FIG. 4: (Color on-line) DOS for type II antiferromagnetic
KNiF3 obtained with a) LDA, b) ASIC1 and c) LDA+U
(U=8 eV and J=1 eV). The valence band top is aligned at
E=0 eV (vertical line). The experimental UPS spectrum from
reference [56] is also presented (thick green line). The relative
binding energy is shifted in order to match the K 3p peak.
ergy part of the band, while the F 2p the lower. For both
the total bandwidth is rather similar and it is about 9-
10 eV. In contrast LDA+U offers a picture where the
Ni-F hybridization spread across the whole bandwidth,
which is now reduced to less than 7 eV.
Experimentally, ultraviolet photoemission spec-
troscopy (UPS) study of the whole KMF3 (M: Mn, Fe,
Co, Ni, Cu, Zn) series [56] gives us a spectrum dominated
by two main peaks: a low K 3p peak and broad band
mainly attributed to F 2p. These two spectroscopical
features are separated by a binding energy of about
10 eV. In addition the 10 eV wide F 2p band has some
fine structure related to various Ni d multiplets. An
analysis based on averaging the multiplet structure [56]
locates the occupied Ni d states at a bounding energy
about 3 eV smaller than that of the F 2p band. In figure
4 we superimpose the experimental UPS spectrum to
our calculated DOS, with the convention of aligning in
each case the sharp K 3p peak.
It is then clear that ASIC provides in general a better
agreement with the UPS data. In particular both the
Ni-3d/F 2p bandwidth and the position of the Fermi en-
ergy (EF) with respect to the K 3p peak are correctly
predicted. This is an improvement over LDA, which de-
scribes well the Ni-3d/F 2p band, but positions the K 3p
states too close to EF. For this reason, when we align the
Method a0 Jth P
d Jex P
LDA 3.951 46.12 (53.1) 1.829 40.4 1.834
PBE 4.052 33.98 (37.0) 1.813 36.48 1.808
BLYP 4.091 31.10 (37.6) 1.821 36.72 1.812
ASIC1/2 3.960 40.83 1.876 36.14 1.878
ASIC1 3.949 36.22 1.907 30.45 1.914
ASIC∗1/2 3.969 43.44 1.876 38.57 1.878
ASIC∗1 3.949 39.80 1.907 33.56 1.914
LDA+U 4.007 12.55 10.47 1.940
TABLE II: Calculated J parameters (in meV) and the
Mülliken magnetic moment for Ni 3d (Pd) in KNiF3. The ex-
perimental values for J and a0 are 8.2±0.6 meV and 4.014Å
respectively while the values in brackets are those from refer-
ence [9]. In the table we report values evaluated at the theo-
retical (Jth and P
d ) and experimental (Jex and P
d ) lattice
constant. ASIC∗1/2 and ASIC
1 are obtained from the LDA
energies evaluated at the ASIC density.
UPS spectrum at the K 3p position, this extends over EF.
Finally in the case of LDA+U , there is a substantial mis-
alignment between the UPS data and our DOS. LDA+U
in fact erroneously pushes part of the Ni d mainfold below
the F 2p DOS, which now forms a rather narrow band.
We now turn our attention to total energy related
quantities. In table II we present the theoretical equi-
librium lattice constant a0 and the Heisenberg exchange
parameter J for all the functionals used. Experimentally
we have J=8.2± 0.6 meV [24]. The values of a0 and J are
calculated for the type II antiferromagnetic ground state,
by constructing a supercell along the (111) direction. Im-
portantly values of J obtained by considering a supercell
along the (100) direction, i.e. by imposing antiferromag-
netic alignment between ferromagnetic (100) planes (type
I antiferromagnet), yield essentially the same result, con-
firming the fact that the interaction is effectively only
extending to nearest neighbors. Furthermore we report
results obtained both at the theoretical equilibrium lat-
tice constant (Jth) and at the experimental one (Jex).
Also in this case local XC functionals largely overesti-
mate J , with errors for Jex going from a factor 8 (LDA)
to a factor 4.5 (GGA-type). ASIC improves these val-
ues, although only marginally, and our best agreement
is found for ASIC1, while ASIC1/2 is substantially iden-
tical to GGA. Interestingly the ASIC1 performance is
rather similar, if not better, to that of meta-GGA func-
tionals [9]. The situation is however worsened when we
consider J parameters obtained at the theoretical lattice
constant. The ASIC-calculated a0 are essentially identi-
cal to those from LDA and about 2% shorter than those
from GGA. Since J depends rather severely on the lattice
parameter we find that at the theoretical lattice constant
GGA-functionals perform actually better than our ASIC.
Finally, also in this case the J ’s obtained by simply us-
ing the LDA energies are larger than those calculated by
including the SIC corrections (see equation 1).
In general the improvement of the J parameter is cor-
related to an higher degree of electron localization, in
particular of the Ni d shell. In table II the magnetic mo-
ment of the Ni d shell Pd, obtained from the Mülliken
population, is reported. This increases systematically
when going from LDA to GGA to ASIC approaching the
atomic value expected from Ni2+.
Our best result is obtained with LDA+U , which re-
turns an exchange of 10.47 meV for the same U and J
that fit the experimental lattice constant. This is some-
how superior performance of LDA+U with respect to
ASIC should not be surprising and it is partially related
to an increased localization. The Ni ions d shell in oc-
tahedral coordination splits into t2g and eg states, which
further split according to Hund’s rule. The t2g states are
all filled, while for the eg only the majority are. By look-
ing at the LDA DOS one can recognize the occupied t↑2g
orbitals (we indicate majority and minority spins respec-
tively with ↑ and ↓) at -3 eV, the e↑g at -2 eV and the t
at about 0 eV, while the empty e↓g are at between 1 and
3 eV above the valence band maximum.
The local Hund’s split can be estimated from the e↑g-
e↓g separation. The ASIC scheme corrects only occupied
states [57], and therefore it enhances the local exchange
by only a downshift of the valence band. From the DOS
of figure 4 it is clear that this is only a small contri-
bution. In contrast the LDA+U scheme also corrects
empty states, effectively pushing upwards in energy the
e↓g band. The net result is that of a much higher degree
of localization of the d shell with a consequent reduction
of the Ni-Ni exchange. This is similar to the situation de-
scribed by the Hartree-Fock method, which however re-
turns exchange parameters considerably smaller than the
experimental value [20, 21, 22, 23]. Interestingly hybrid
functionals [7] have the right mixture of non-local ex-
change and electron correlation and produce J ’s in close
agreement with the experiments.
We further investigate the magnetic interaction by
evaluating J as a function of the lattice constant. Ex-
perimentally this can be achieved by replacing K with
Rb and Tl, and indeed de Jongh and Block [58] early
suggested a d−α power law with α = 12± 2. Our calcu-
lated J as a function of the lattice constant d for LDA,
GGA, ASIC1 and LDA+U (U=8 eV and J=1 eV) are
presented in figure ??. For all the four functionals in-
vestigated J varies as a power law, although the calcu-
lated exponents are rather different: 8.6 for LDA, 9.1
for GGA, 11.3 for ASIC1 and 14.4 for LDA+U . This
further confirms the strong underestimation of the ex-
change constants from local functionals. Clearly the rel-
ative difference between the J obtained with different
functionals becomes less pronounced for small d, where
the hybridization increases and local functionals perform
better. Note that only ASIC1 is compatible with the ex-
perimental exponent of 12 ± 2, being the one evaluated
from LDA+U too large. Importantly we do not expect
to extrapolate the LDA+U value at any distance, since
FIG. 5: J as a function of the lattice constant for LDA, GGA,
ASIC1 and LDA+U (U=8 eV and J=1 eV). The symbols are
our calculate value while the solid lines represent the best
power-law fit.
the screening of the parameters U and J changes with
the lattice constant.
In conclusion for the critical case of KNiF3 the ASIC
method appears to improve the LDA results. This is es-
sentially due to the better degree of localization achieved
by the ASIC as compared with standard local function-
als. However, while the improvement over the bandstruc-
ture is substantial, it is only marginal for energy-related
quantities. The main contribution to the total energy in
the ASIC scheme comes from the LDA functional, which
is now evaluated at the ASIC density. This is not suf-
ficient for improving the exchange parameters, which in
contrast need at least a portion of non-local exchange.
C. Transition metal monoxides
Another important test for the ASIC method is that of
transition metal monoxides. These have been extensively
studied both experimentally and theoretically and they
are the prototypical materials for which the LDA appears
completely inadequate. In this work we consider MnO
and NiO, which have respectively half-filled and partially-
filled 3d shells. They both crystallize in the rock-salt
structure and in the ground state they are both type-
II antiferromagnetic insulators. The Néel’s temperatures
are 116 K and 525 K respectively for MnO and NiO. In all
our calculations we consider double-ζ polarised basis for
the s and p shell of Ni, Mn and O, double-ζ for the Ni and
Mn 3d orbitals, and single-ζ for the empty 3d of O. We
sample 6×6×6 k-points in the full Brillouin zone of both
the cubic and rhombohedral cell describing respectively
type I and type II antiferromagnetism. Finally the real-
space mesh cutoff is 500 Ry.
The calculated band structures obtained from LDA,
ASIC1/2 and ASIC1 are shown in figures 6 and 7 for MnO
and NiO respectively. These have been already discussed
extensively in the context of the ASIC method [27, 28]
and here we report only the main features. For both the
materials LDA already predicts an insulating behavior,
although the calculated gaps are rather small and the
nature of the gaps is not what experimentally found. In
both cases the valence band top has an almost pure d
component, which suggests these materials to be small
gap Mott-Hubbard insulators. The ASIC downshifts the
occupied d bands which now hybridize with the O-p man-
ifold. The result is a systematic increase of the band-gap
which is more pronounced as the parameter α goes from
1/2 to 1. Importantly, as noted already before [28], the
experimental band-gap is obtained for α ∼ 1/2.
FIG. 6: Calculated band structure for the type II anti-
ferromagnetic MnO obtained from a) LDA, b) ASIC1/2 and
c) ASIC1. The valence band top is aligned at 0 eV (horizontal
line).
We then moved to calculating the exchange parame-
ters. In this case we extend the Heisenberg model to
second nearest neighbors, by introducing the first (J1)
and second (J2) neighbor exchange parameters. These
are evaluated from total energy calculations for a ferro-
magnetic and both type II and type I antiferromagnetic
alignments. Our calculated results, together with a few
selected data available from the literature are presented
in table III.
Let us first focus our attention to MnO. In this case
both the J ’s are rather small and positive (antiferro-
magnetic coupling is favorite), in agreement with the
Goodenough-Kanamori rules [59] and the rather low Néel
temperature. Direct experimental measurements of the
exchange parameters are not available and the com-
monly accepted values are those obtained by fitting the
magnetic susceptibility with semi-empirical methods [60].
Importantly this fit gives almost identical first and second
nearest neighbour exchange constants. In contrast all the
exchange functionals we have investigated offer a picture
where J2 is always approximately twice as large as J1.
FIG. 7: Calculated band structure for the type II anti-
ferromagnetic NiO obtained from a) LDA, b) ASIC1/2 and
c) ASIC1. The valence band top is aligned at 0 eV (horizon-
tal line).
Method MnO NiO
J1 J2 Pd J1 J2 Pd
LDA 1.0 2.5 4.49 (4.38) 13.0 -94.4 1.41 (1.41)
PBE 1.5 2.5 4.55 (4.57) 7.0 -86.8 1.50 (1.59)
ASIC1/2 1.15 2.44 4.72 (4.77) 6.5 -67.3 1.72 (1.77)
ASIC1 0.65 1.81 4.84 (4.86) 3.8 -41.8 1.83 (1.84)
ASIC∗1/2 1.27 2.65 4.72 (4.77) 7.1 -74.6 1.72 (1.77)
ASIC∗1 0.69 2.03 4.84 (4.86) 4.4 -47.9 1.83 (1.84)
SE1a 0.86 0.95
HFb 0.22 0.36
B3LYPc 0.81 1.71
PBE0b 0.96 1.14
B3LYPd 2.4 -26.7
HFd 0.8 -4.6
SIC-LDAe 2.3 -12
Expt.f 1.4 -19.8
Expt.g 1.4 -17.0
TABLE III: Calculated J1 and J2 in meV for MnO and NiO.
Pd is the magnetic moment of the d shell calculated from the
type II antiferromagnetic phase. Values in bracket are for Pd
evaluated from the ferromagnetic ground state. ASIC∗1/2 and
ASIC∗1 are obtained from the LDA energies evaluated at the
ASIC density. a) Ref. [60], b) Ref. [61], c) Ref. [62], d) Ref.
[11], e) Ref. [25], f) Ref. [64], g) Ref. [65]
This gives us a reasonably accurate value of J1 for LDA
and GGA, but J2 is overestimated by approximately a
factor 2, in agreement with previous calculations [10].
ASIC systematically improves the LDA/GGA descrip-
tion, by reducing both J1 and J2. This is related to the
enhanced localization of the Mn d electrons achieved by
the ASIC, as it can be seen by comparing the Mn d mag-
netic moments (Pd) calculated for different functionals
(see table III). Thus ASIC1, which provides the largest
magnetic moment, gives also J ’s in closer agreement with
the experimental values, while ASIC1/2 is not very dif-
ferent from LDA.
Importantly for half-filling, as in MnO, the ASIC
scheme for occupied states is fundamentally analogous
to the LDA+U method, with the advantage that the U
parameter does not need to be evaluated. Finally, at
variance with KNiF3, it does not seem that a portion of
exact exchange is strictly necessary in this case. Hartree-
Fock [61] results in a dramatic underestimation of the J
parameters, while B3LYP [62] is essentially very similar
to ASIC1. Curiously the best results available in the lit-
erature [61] are obtained with the PBE0 functional [63],
which combines 25% of exact-exchange with GGA.
The situation for NiO is rather different. The ex-
perimentally available data [64, 65] show antiferromag-
netic nearest neighbour and ferromagnetic second near-
est neighbour exchange parameters. The magnitude is
also rather different with |J2| > 10 |J1|. Standard lo-
cal functionals (LDA and GGA) fail badly and overes-
timate both the J ’s by more than a factor 5. ASIC
in general reduces the exchange constants and drasti-
cally improves the agreement between theory and exper-
iments. In particular ASIC1 gives exchange parameters
only about twice as large as those measured experimen-
tally.
A better understanding can be obtained by looking at
the orbital-resolved DOS for the Ni d and the O p orbitals
(figure 8) as calculated from LDA and ASIC. There are
two main differences between the LDA and the ASIC
results. First there is an increase of the fundamental
band-gap from 0.54 eV for LDA to 3.86 eV for ASIC1/2 to
6.5 eV for ASIC1. Secondly there is change in the relative
energy positioning of the Ni d and O p contributions to
the valence band. In LDA the top of the valence band is
Ni d in nature, with the O p dominated part of the DOS
lying between 4 eV and 8 eV from the valence band top.
ASIC corrects this feature and for ASIC1/2 the O p and
Ni d states are well mixed across the whole bandwidth. A
further increase of the ASIC corrections (α = 1) leads to
a further downshift of the Ni d band, whose contribution
becomes largely suppressed close to the valence band-top.
Thus, increasing the portion of ASIC pushes NiO further
into the charge transfer regime.
Interestingly, although ASIC1/2 gives the best band-
structure, the J ’s obtained with ASIC1 are in better
agreement with the experiments. This is somehow sim-
ilar to what observed when hybrid functionals are put
to the test. Moreira et al. demonstrated [11] that J ’s
in close agreement with experiments can be obtained by
using 35% Hartree-Fock exchange in LDA. Moreover, in
close analogy to the ASIC behaviour, as the fraction of
exact exchange increases from LDA to purely Hartree-
Fock, the exchange constants decrease while the band-
gap gets larger. However, while the best J ’s are obtained
with 35% exchange, a gap close to the experimental one
FIG. 8: Calculated orbital resolved DOS for type II anti-
ferromagnetic NiO obtained with a) LDA, b) ASIC1/2 and c)
ASIC1. The valence band top is aligned at 0 eV.
is obtained with B3LYP, which in turns overestimates the
J ’s. This remarks the subtile interplay between exchange
and correlations in describing the magnetic interaction
of this complex material. Finally, it is worth remarking
that a fully self-consistent SIC [25] seems to overcorrect
the J ’s, while still presenting the erroneous separation
between the Ni d and O p states.
IV. CONCLUSIONS
In conclusions the approximated expression for the
ASIC total energy is put to the test of calculating ex-
change parameters for a variety of materials, where local
and gradient-corrected XC functionals fail rather badly.
This has produced mixed results. On the one hand, the
general bandstructure and in particular the valence band,
is considerably improved and resembles closely data from
photo-emission. On the other hand, the exchange con-
stants are close to experiments only for the case when the
magnetism originates from half-filled shells. For other
fillings, as in the case of NiO or KNiF3 the ASIC im-
provement over LDA is less satisfactory, suggesting that
a much better XC functional, incorporating a portion at
least of exact exchange, is needed. Importantly ASIC
seems to be affected by the same pathology of hybrid
functional, i.e. the amount of ASIC needed for correct-
ing the J is different from that needed for obtaining a
good bandstructure.
V. ACKNOWLEDGEMENTS
This work is supported by Science Foundation of Ire-
land under the grant SFI05/RFP/PHY0062. Computa-
tional resources have been provided by the HEA IITAC
project managed by the Trinity Center for High Perfor-
mance Computing and by ICHEC.
[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864
(1964).
[2] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[3] W. Kock and M. C. Holthausen, A Chemist’s Guide
to Density Functional Theory, (Wiley-VCH, Weinheim,
2000).
[4] I. Turek, J. Kudronovsky, V. Drchals and P. Bruno, Phil.
Mag. 86, 1713 (2006).
[5] E. Ruiz, S. Alvarez, J. Cano, and V. Polo, J. Chem. Phys.
123, 164110 (2005).
[6] E. Ruiz, J. Cano, S. Alvarez, and P. Alemany, J. Comput.
Chem. 20, 1391 (1999).
[7] R. L. Martins and F. Illas, Phys. Rev. Lett. 79, 1539
(1997).
[8] F. Illas and R. L. Martins, J. Chem. Phys. 108, 2519
(1998).
[9] I. Ciofini, F. Illas and C. Adamo, J. Chem. Phys. 120,
3811 (2004).
[10] J. E. Pask, D. J. Singh, I. I. Mazin, C. S. Hellberg and
J. Kortus, Phys. Rev. B 64, 024403 (2001).
[11] I. de P. R. Moreira, F. Illas, and R. Martin, Phys. Rev.
B 65, 155102 (2002).
[12] E. Ruiz, A. Rodriguez-Fortea, J. Cano, S. Alvarez and P.
Alemany, J. Comput. Chem. 24, 982 (2003).
[13] J. Cabrero, N. Ben Amor, C. de Graaf, F. Illas and R.
Caballol, J. Phys. Chem. A 104, 9983 (2000).
[14] A. Bencini, F. Totti, C. A. Daul, K. Doclo, P. Fantucci,
and V. Barone, Inorg. Chem. 36, 5022 (1997).
[15] L. Noodleman, J. Chem. Phys. 74, 5737 (1981).
[16] K. Yosida, Theory of Magnetism, Vol. 122 Springer Series
in Solid State Sciences (Springer, Heidelberg, 1992).
[17] R. G. Parr and W.Yang, Density functional theory of
atoms and molecules, (Oxford University Press, 1989).
[18] J. Tao, J.P. Perdew, V.N. Staroverov and G.E. Scuseria,
Phys. Rev. Lett. 91, 146401 (2003).
[19] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048
(1981).
[20] J. M. Ricart, R. Dovesi, C. Roetti and V. R. Saunders,
Phys. Rev. B 52, 2381 (1995).
[21] R. Dovesi, J. M. Ricart, V. R. Saunders and R. Orlando,
J. Phys.: Condens. Matter 7, 7997 (1995).
[22] I. de P. R. Moreira and F. Illas, Phys. Rev. B 55, 4129
(1997).
[23] R. Caballol, O. Castell, F. Illas, I. de P. R. Moreira and
J. P. Malrieu, J. Phys. Chem. A 101, 7860 (1997).
[24] M. E. Lines, Phys. Rev. 164, 736 (1967).
[25] D. Ködderitzsch, W. Hergert, W.M. Temmerman,
Z. Szotek, A. Ernst and H. Winter, Phys. Rev. B 66,
064434 (2002).
[26] D. Vogel, P. Kruger and J. Pollamann, Phys. Rev. B 54,
5495 (1996).
[27] A. Filippetti and N. A. Spaldin, Phys. Rev. B 67, 2435
(2003).
[28] C. D. Pemmaraju, T. Archer, D. Sanchez-Portal and S.
Sanvito, Phys. Rev. B 75, 045101 (2007).
[29] A. Filippetti and V. Fiorentini, Phys. Rev. Lett. 95,
086405 (2005).
[30] D. Vogel, P. Krüger and J. Pollmann, Phys. Rev. B 58,
3865 (1998).
[31] A. Filippetti and V. Fiorentini, Rev. B 73, 035128 (2005).
[32] A. Filippetti and N.A. Spaldin, Phys. Rev. B 68, 045111
(2003).
[33] B.B. Van Aken, T.T.M.A. Palstra, A. Filippetti and
N.A. Spaldin, Nature Mat. 3, 164 (2003).
[34] P. Delugas, V. Fiorentini and A. Filippetti, Phys. Rev.
B 71, 134302 (2005).
[35] A. Filippetti, N.A. Spaldin and S. Sanvito, Chem. Phys.
309, 59 (2005).
[36] A. Filippetti, N.A. Spaldin and S. Sanvito, J. Magn.
Magn. Mater. 290, 1391 (2005).
[37] C. Toher, A. Filippetti, S. Sanvito and K. Burke, Phys.
Rev. Lett. 95, 146402 (2005).
[38] C. Toher and S. Sanvito, cond-mat/0611617
[39] Note that here we use “LDA” also to refer to the local
spin density approximation (LSDA), i.e. to the spin un-
restricted version of the theory.
[40] J.G. Harrison, R.A. Heaton and C.C. Lin, J. Phys. B:
At. Mol. Phys. 16, 2079 (1983).
[41] M.R. Pederson, R.A. Heaton and C.C. Lin, J. Chem.
Phys. 80, 1972 (1984).
[42] M.R. Pederson, R.A. Heaton and C.C. Lin, J. Chem.
Phys. 82, 2688 (1985).
[43] A. Svane and O. Gunnarsson, Phys. Rev. Lett. 65, 1148
(1990).
[44] Z. Szotek, W.M. Temmerman and H. Winter, Phys. Rev.
B 47, 4029 (1993).
[45] J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Jun-
quera, P. Ordejon, and D. Sanchez-Portal, J. Phys.: Con-
dens. Matter 14, 2745 (2002).
[46] D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566
(1980).
[47] A. D. Becke, Phys. Rev. A 38, 3098 (1988).
[48] C. Lee, W. Yang and R. G. Parr, Phys. Rev. B 37, 785
(1988).
[49] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev.
Lett. 77, 3865 (1996).
[50] C. Adamo, V. Barone, A. Bencini, R. Broer, M. Filatov,
N.M. Harrison, F. Illas, J.P. Malrieu, and I. de P.R. Mor-
eira, J. Chem. Phys. 124, 107101 (2006).
[51] E. Ruiz, S. Alvarez, J. Cano, and V. Polo, J. Chem. Phys.
http://arxiv.org/abs/cond-mat/0611617
124, 107102 (2006).
[52] L. J. de Jongh and R. Miedema, Adv. Phys. 23, 1 (1974).
[53] V. I. Anisomov, M. A. Korotin, J. A. Zaanen, and O. K.
Anderson, Phys. Rev. Lett. 68, 345 (1992).
[54] M. T. Czyzyk and G. A. Sawatzsky, Phys. Rev. B 40,
14211 (1994).
[55] M. Wierzbowska, D. Sánchez-Portal and S. Sanvito,
Phys. Rev. B 70, 235209 (2004).
[56] H. Onuki, F. Sugawara, M. Hirano and Y. Yamaguchi, J.
Phys. Soc. Jpn. 49, 2314 (1980).
[57] The ASIC method is not completely self-interaction free
and therefore some spurious corrections to the unoccu-
plied KS states are present.
[58] L.J. de Jongh and R. Block, Physica 79B, 568 (1975).
[59] J. Goodenough, Magnetism and the Chemical Bond,
(John Wiley and Sons, New York, 1963).
[60] M. E Lines and E. D. Jones, Phys. Rev. 139, A1313
(1965).
[61] C. Franchini, V. Bayer, R. Podloucky, J. Paier and
G. Kresse, Phys. Reb. B 72, 045132 (2005).
[62] X. Fenf, Phys. Rev. B 69, 155107 (2004).
[63] M. Ernzerhof and G.E. Scuseria, J. chem. Phys. 110,
5029 (1999).
[64] M. T. Hutching and E. J. Samuelsen, Phys. Rev. B 6,
3447 (1972).
[65] R. Shanker and R. A. Singh, Phys. Rev. B 7, 5000 (1973).
[66] S. Hüfner and T. Riesterer, Phys. Rev. B 33, 7267 (1986).
	Introduction
	The atomic SIC method
	Results
	H-He chain
	Ionic antiferromagnets: KNiF3
	Transition metal monoxides
	Conclusions
	Acknowledgements
	References
