Ab initio wavefunction based methods for excited states in solids:
  correlation corrections to the band structure of ionic oxides
Ab initio wavefunction based methods for excited states in solids: correlation
corrections to the band structure of ionic oxides
L. Hozoi, U. Birkenheuer, and P. Fulde
Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Str. 38, 01187 Dresden, Germany
A. Mitrushchenkov and H. Stoll
Universität Stuttgart, Pfaffenwaldring 57, 70550 Stuttgart, Germany
(Dated: November 5, 2018)
Ab initio wavefunction based methods are applied to the study of electron correlation effects on
the band structure of oxide systems. We choose MgO as a prototype closed-shell ionic oxide. Our
analysis is based on a local Hamiltonian approach and performed on finite fragments cut from the
infinite solid. Localized Wannier functions and embedding potentials are obtained from prior peri-
odic Hartree-Fock (HF) calculations. We investigate the role of various electron correlation effects
in reducing the HF band gap and modifying the band widths. On-site and nearest-neighbor charge
relaxation as well as long-range polarization effects are calculated. Whereas correlation effects are
essential for computing accurate band gaps, we found that they produce smaller changes on the
HF band widths, at least for this material. Surprisingly, a broadening effect is obtained for the O
2p valence bands. The ab initio data are in good agreement with the energy gap and band width
derived from thermoreflectance and x-ray photoemission experiments. The results show that the
wavefunction based approach applied here allows for well controlled approximations and a transpar-
ent identification of the microscopic processes which determine the electronic band structure.
I. INTRODUCTION
The proper treatment of electron correlation effects in
molecules and solids stands at the heart of modern elec-
tronic structure theory1,2,3,4,5. In the study of molecular
systems at least, wavefunction based quantum chemistry
is known to provide a rigorous theoretical framework for
addressing the electron correlation problem. The stan-
dard quantum chemical methods4,5 make possible the
construction of approximate wavefunctions at levels of
increasing sophistication and accuracy and offer thus a
systematic route to converged results. Advanced wave-
function based calculations can be routinely performed
nowadays for small and medium size molecules. However,
algorithms able to treat electron correlation effects in pe-
riodic systems are still at their infancy. The simplest cor-
relation method is based on the Møller-Plesset equations
and second-order Møller-Plesset (MP2) perturbational
schemes for solids were actually implemented by several
groups6,7,8,9,10,11,12,13,14. As an extension of MP2, ab ini-
tio many-body Green’s function techniques were devel-
oped too15,16,17,18,19. In addition, investigations based on
coupled-cluster (CC) theory were initiated13,20,21,22,23.
The main point when applying quantum chemical
methods to correlation calculations is to make use of the
local character of the correlation hole which is surround-
ing an electron. The latter optimizes the Coulomb re-
pulsion between electrons and its accurate description is
the essence of the correlation treatment. Starting point
is a Hartree-Fock (HF) calculation on top of which the
correlation calculations are implemented. For insulators
and semiconductors, the development of efficient compu-
tational tools is facilitated by the use of optimally local-
ized Wannier functions. Several orbital localization pro-
cedures were proposed in the context of periodic HF cal-
culations. They normally rest on the a posteriori trans-
formation of the crystal Bloch orbitals, see for example
Refs. 6,14,20,21,22,24,25. Alternative approaches were
elaborated by Shukla et al.26,27,28,29 and in the Toulouse
group11,30, where the self-consistent-field (SCF) equa-
tions are solved directly in the Wannier representation.
Local correlation methods such as the incremental
scheme of Stoll31,32 and the local Hamiltonian formal-
ism of Horsch et al.33 and Gräfenstein et al.34,35,36
were applied before to the rigorous determination of
both ground-state13,26,27,28 and excited-state18,37,38,39,40
properties of infinite systems. Materials under inves-
tigation were crystalline LiH, LiF, and LiCl18,26,27,28,
diamond38,40, silicon38, beryllium39, carbon, and boron-
nitrogen chains13, and trans-polyacetylene37. In the
present paper we extend these studies to the case of ionic
oxide compounds and choose MgO as a prototype and,
at the same time, relatively simple insulating oxide. We
compute correlation induced corrections to the valence
and conduction HF energy bands and find that our final
estimate for the fundamental gap is in rather good agree-
ment with the experimental data. Our analysis provides
also a clear picture of the major correlation effects that
are responsible for the reduction of the HF band gap.
The present study on magnesium oxide should open the
way to similar calculations for more complex materials
such as the 3d transition-metal oxides.
II. COMPUTATIONAL APPROACH
Periodic HF calculations were carried out for crys-
talline MgO with the crystal package41. Wannier
functions associated with the HF valence and low-lying
conduction bands were determined with the Wannier-
http://arxiv.org/abs/0704.1906v1
Boys localization module24 of the same program. For
conduction-band states, this may be a somewhat tedious
task. However, localized Wannier orbitals can also be de-
rived in more difficult cases with band disentanglement
techniques25.
For computing the correlation induced corrections to
the valence and conduction bands we adopt the same
quasiparticle picture and local Hamiltonian formalism as
employed in Refs. 33,34,35,36,37,38,39,40. The physical
reasoning is as follows. When an electron (or hole) is
added to the N -particle system, its surroundings relax
and polarize due to the additional charge. This response
of the system lowers the energy it takes to add the ex-
tra electron (or hole). The particle plus the modified
surroundings move together through the system in the
form of a Bloch wave and define a quasiparticle. At the
same time, some of the correlation contributions which
are present in the ground-state of the N -particle system
are no longer operative in the (N+1) [or (N−1)] -particle
system, because some excitations are now blocked. This
effect is called loss of ground-state correlation and must
be accounted for too33,38. A basic feature of our approach
is that the correlation treatment is performed on a finite
cluster C cut from the extended periodic system. With
the cluster-in-solid embedding technique elaborated by
Birkenheuer et al.38,42,43, the cluster C is divided into an
active region CA that supports the local occupied and
virtual orbitals entering the post-HF calculations and a
spatial “buffer” domain CB including a number of atomic
sites (and basis functions) whose role is to provide a good
representation for the tails of the Wannier-like orbitals of
the active region CA. Thus C = CA + CB. All Wannier
orbitals centered in CB (and also the environment) will
be held frozen in the correlation calculations.
The orbital set associated with the finite-size clus-
ter is obtained from the data supplied by the Wannier-
Boys localization module of crystal24. The interface
program44 written for this purpose yields in addition an
embedding potential corresponding to the frozen environ-
ment, i.e., the surrounding HF electron sea. The detailed
procedure for constructing the embedded cluster and the
associated orbital set is described elsewhere42,43. Here
we merely give an outline of this. The core, valence, and
low-lying conduction-band cluster orbitals are generated
by projecting the crystal Wannier functions onto the set
of atomic basis functions attached to the cluster region
|w′n(R)〉 =
β,α∈C
|β〉S−1
〈α|wn(R)〉 , (1)
where α and β are Gaussian basis functions centered
within the region C, S−1
represents the inverse overlap
matrix for the basis set attached to C, |wn(R)〉 is a Wan-
nier orbital with index n and centered in the unit cell
corresponding to the lattice vector R, and |w′n(R)〉 is
its projected counterpart expressed exclusively in terms
of basis functions centered in C. The |w′n(R)〉 func-
tions are neither normalized nor are they orthogonal to
each other because of the projection procedure mentioned
above. Therefore they are group-wise orthonormalized in
the following order: active core, active occupied, active
low-lying conduction-band orbitals, buffer core, buffer oc-
cupied, buffer low-lying conduction-band orbitals. This
way the contamination of the most important types of or-
bitals is minimized. Group-wise orthonormalization ac-
tually means to use Schmidt orthogonalization for the
inter-group orthogonalization and Löwdin orthonormal-
ization inside the groups. This set of orthonormal or-
bitals will be denoted by |w̃′n(R)〉. For the construction
of the variational space to be used in the subsequent cor-
relation calculations we follow the prescription suggested
by Pulay and Saebø45. So-called virtual projected atomic
orbitals (PAO’s) are generated from the Gaussian basis
functions associated with the active region CA by pro-
jecting out the occupied and the low-lying conduction-
band orbitals |w̃′n(R)〉 via a Schmidt orthogonalization
scheme. Thereafter, the PAO’s are Löwdin orthonormal-
ized among themselves to facilitate their subsequent use
in the correlation calculations.
By making the buffer region CB sufficiently large the
original Wannier orbitals of CA are well represented. For
example, for the Wannier orbitals corresponding to the
lowest four conduction bands of MgO and centered at
a Mg site, the six ligands surrounding the given Mg ion
should be included in CA because large contributions arise
not only from the Mg 3s and 3p basis functions but also
from the nearest-neighbor (NN) oxygen 2s,3s and 2p,3p
components. Farther neighbors need not be included in
the central region CA but may be put into the buffer
zone CB because the weight of the longer-range tails is
small. Whereas a high-quality description can be easily
achieved for the Wannier orbitals centered in the active
region, the representation of the Wannier functions in the
buffer zone is less accurate. The impact of this deficiency
on the correlation calculations, however, is completely
compensated by an appropriate choice of the embedding
potential. Actually, the Gaussian orbital representation
V embαβ of the embedding potential is constructed by
V embαβ = F
− F [P C ]αβ , α, β ∈ C (2)
where F crys is the self-consistent Fock operator from the
periodic HF calculation and F [P C ] is the Fock operator
associated with the density operator
P C = 2
|w̃′ν〉〈w̃′ν | (3)
arising from all occupied orbitals |w̃′ν〉 which enter the
subsequent correlation calculations explicitly. This way,
the correlation calculations are effectively performed in
an infinite frozen HF environment.
The data concerning the occupied and virtual orbitals
of the cluster is transferred via the crystal-molpro
interface44 to the quantum chemistry program mol-
46. The same holds for the matrix representation
of the self-consitent Fock operator of the periodic
host system. The embedding potential itself is con-
structed according to Eqs. (2) and (3) using the MA-
TROP module of the molpro program package.
Local electron removal and electron addition one-
particle configurations can be defined in terms of the set
of occupied and virtual orbitals localized within the spa-
tial domain CA:
|ΦN−1
〉 = cRpσ|Φ〉 and |ΦN+1Rqσ 〉 = c
|Φ〉 , (4)
where cRpσ and c
are annihilation and creation op-
erators for the valence and conduction band σ-spin or-
bitals |w̃′p(R)〉 and |w̃′q(R)〉, respectively, and |Φ〉 is the
single-determinant ground-state wavefunction of the N -
electron system. For clusters which are large enough,
the Hartree-Fock valence and conduction energy bands
of the periodic crystal can be recovered by diagonalizing
k-dependent matrices of the following form:
HHFnn′(k) =
eikR〈ΦN∓1
0nσ |H − E
0 |ΦN∓1Rn′σ〉 . (5)
Here 0 stands for the reference unit cell and EHF0 denotes
the ground-state HF energy of the neutral N -electron
system. The diagonal terms R=0, n′=n in this expres-
sion are directly related to the on-site excitation energies
within the Koopmans approximation4,5, i.e., ionization
potentials
IPHFpp (0) = 〈ΦN−10pσ |H |Φ
0pσ 〉 − E
0 = −ǫHF0p > 0 (6)
and electron affinities
EAHFqq (0) = E
0 − 〈ΦN+10qσ |H |Φ
0qσ 〉 = −ǫ
0q , (7)
with the latter being negative in the case of MgO. The
off-diagonal terms are the hopping matrix elements in a
tight-binding representation,
tHFnn′(R) = 〈ΦN∓10nσ |H − EHF0 |Φ
〉 . (8)
Starting from the Hartree-Fock energy bands we in-
clude next the effects of electron correlations. Within
the quasiparticle approximation we may introduce for the
(N+1)-particle system a wavefunction of the following
form:
|ΨN+1
〉 = exp (S) c†
|Ψ〉 , (9)
where |Ψ〉 is the ground-state of the N electron system.
The operator S does not commute with c
and de-
scribes the correlation hole. For a more detailed discus-
sion we refer to Ref. 1. Alternatively, we may write
|ΨN+1
〉 = Ω c†
|Φ〉, where |Φ〉 is the SCF ground-state
and Ω acts like a wave- or Moeller operator3. Therefore,
when correlations are taken into account, the analogue
of (5) becomes
Hnn′(k) =
eikR〈ΨN∓1
0nσ |H − E0|Ψ
〉 , (10)
where E0 is the energy of the N -particle correlated
ground-state |Ψ〉 (see also Refs. 34,35,36,37,38). The cor-
relation hole around an added electron (or hole) consists
of a short-range and a long-range part. The short-range
part originates from intra-atomic and short-range inter-
atomic relaxation and polarization effects. We construct
the short-range part of the correlation hole by separate
orbital optimizations. Thereby the Wannier orbital to
which the extra electron (hole) is attached is kept frozen
and the changes within a finite region around it are deter-
mined by performing an additional SCF calculation38,39.
This is sufficient if electron correlations are weak or mod-
erate. If they are strong we would have to take into
account the fact that changes in the nearby surround-
ings are decreased when correlation effects among the
electrons in that neighborhood are accounted for. The
long-range part of the correlation hole consists of long
range polarization of the environment. The effect of long-
range polarization on the diagonal Hamiltonian matrix
elements is estimated in this paper by applying the ap-
proximation of a dielectric continuum. With the known
dielectric constant of MgO, ǫ0=9.7, we can directly de-
termine the polarization energy of a charge ±e outside a
sphere of radiusR. In addition, differential correlation ef-
fects related to the existence of a different number of elec-
trons in the system’s ground-state and in the (N±1) ex-
cited states are investigated by subsequent configuration-
interaction (CI) calculations.
III. CORRELATION CORRECTIONS TO THE
BAND STRUCTURE OF MGO
Magnesium oxide is a prototype closed-shell ionic ma-
terial that crystallizes under normal conditions in the
rocksalt structure. It is extensively used in materials sci-
ence as substrate for the epitaxial growth of films of other
compounds. In recent years, it has attracted renewed in-
FIG. 1: Hartree-Fock band structure of bulk MgO. The core
O 1s and Mg 1s− 2p bands are not shown in the figure.
FIG. 2: Plot of the Mg 3s-like conduction-band Wannier or-
bital after projection onto a [Mg19O14] cluster. There is sub-
stantial weight at the NN oxygen sites.
terest because of its use as tunnel barrier in magnetic tun-
nel junctions47. A different area where MgO attracted at-
tention is optoelectronics. It turns out that when alloyed
with ZnO, depending on the precise chemical composi-
tion, the band gap of the system can be tuned along an
interval ranging from 3.3 to 7.8 eV48. The latter number
represents the fundamental gap of MgO49.
The valence bands of MgO have oxygen 2p character,
whereas the lowest conduction bands are mainly related
to the Mg 3s and 3p orbitals, with some admixture from
the O 2s,3s and 2p,3p functions. The Hartree-Fock en-
ergy bands are shown in Fig.1. We employed for our
calculations the lattice constant reported by Sasaki et
al., a=4.217 Å50, and Gaussian-type basis sets from the
standard crystal library. Basis sets of triple-zeta qual-
ity, Mg 8-51151, were used for the cations and triple-zeta
basis sets suplemented with polarization functions, O 8-
411∗51, were applied for the more polarizable O ions. At
the Hartree-Fock level and with this choice of the basis
functions, the fundamental gap of the system is 16.20
eV, i.e., 8.4 eV larger than the experimental value. As
illustrated in Fig.1, there is a separation of about 16 eV
between the O 2s and O 2p bands. There is also a clear
separation between the low-lying Mg 3s, 3p complex and
the other conduction bands. Density-functional calcula-
tions within the local density approximation and using
the same Gaussian basis sets as for the HF calculations
give a gap of 5.0 eV between the valence and conduction
bands, nearly 3 eV lower than observed in experiment.
Projected Wannier orbitals associated with the Mg 3s,
3p conduction-band complex are plotted in Fig.2 and
Fig.3. We note that the Wannier-Boys localization mod-
ule of the crystal program yields a set of Mg sp3 hybrids
for the lowest four conduction-band states. In order to
arrive at a set of s and x, y, z -like functions, we applied
to those projected hybrids a Pipek-Mezey localization
procedure52, as implemented in the molpro package46.
In this localization scheme, the number of atomic or-
bitals contributing to a given molecular-like composite
is minimized. Whereas the Mg 3s and 3p Wannier func-
tions have substantial weight at the NN ligand sites, the
valence-band O 2p (and 2s) Wannier orbitals are much
more compact and contributions from neighboring ions
are not visible in plots like those shown in Figs. 2 and
3. In the calculations reported here the norm of the pro-
jected active orbitals is typically larger than 99.0% and
never below 98.5% of the original crystal Wannier func-
tions.
A. Correlation induced corrections to the diagonal
matrix elements
Relaxation and polarization effects in the immedi-
ate neighborhood of an oxygen hole were computed by
separate restricted open-shell HF (ROHF) calculations
on [O39Mg30] clusters. The multiconfiguration MCSCF
module of the molpro package was employed for this
purpose. The active region CA of the [O39Mg30] clus-
ter contains a central 2p5 (or 2s1) O− ion and four of the
twelve nearest O2− ligands. These four oxygen neighbors
are all chosen to be in the same plane and we denote them
as O1xy,...O
xy (see Fig.4). In addition, we include in the
cluster C the NN cations and all the O’s in the next co-
ordination shell of each of the Oixy sites, as sketched in
Fig.4. These additional Mg and O neighbors represent
the so-called buffer region CB and ensure an accurate de-
scription of the tails of the orbitals centered in the active
region CA.
When performing the ROHF calculations, the oxygen
hole orbital is kept frozen39. We also freeze in our calcu-
lations the core-like 1s, 2s, and 2p shells of all Mg ions.
Orbital relaxation effects are listed in Table I for both 2s
FIG. 3: Mg p-like conduction-band Wannier orbital after pro-
jection onto a [Mg19O14] cluster.
and 2p oxygen holes. The on-site relaxation effect is quite
large, more than 2 eV. There is also a substantial relax-
ation/polarization effect associated with the first oxygen
neighbors. For the four Oixy ligands included in the ac-
tive region CA of the cluster, this effect amounts to about
0.41 eV in the presence of an O 2s hole and to 0.40 eV
in the presence of a 2p hole. In Table I, we multiplied
these numbers by three because there is a total of twelve
oxygens in that coordination shell. We obtain thus a
good estimate for the relaxation and polarization effects
up to the nearest O neighbors. The fact that these relax-
ation/polarization effects are additive was checked by ex-
tra calculations with smaller, double-zeta basis sets51 on
two different clusters: a [O39Mg30] cluster including only
the four Oixy sites in the active region and a [O55Mg38]
cluster where all twelve first oxygen neighbors were al-
lowed to polarize.
The four lowest-energy (N+1) conduction-band states
imply Mg 3s1 and 3p1 electron configurations. ROHF
calculations were performed for such configurations on a
[Mg19O38] cluster with a [MgO6] kernel as active region,
CA. Beyond the [MgO6] kernel, this cluster incorporates
again all Mg and O ions in the first two coordination
shells of the active ligands. The on-site relaxation effects
associated with the addition of an electron in a localized
Mg 3s or 3p Wannier orbital are vanishingly small. The
relaxation effects at the adjacent O sites induce energy
shifts of 0.80–0.85 eV, see Table I. In these calculations
the open-shell active orbitals (Mg 3s1 or 3p1) were again
kept frozen39. We note that the energetic effect is nearly
the same for the 3s1 and 3p1 conduction-band states. At
the scale of Fig.1 at least, it induces an uniform down-
FIG. 4: Sketch of the [O39Mg30] cluster employed for the
calculation of short-range relaxation effects on the O valence-
band states. The shortest line segments are Mg–O “bonds”.
The so-called active region CA includes a central O ion and
four nearest oxygen neighbors Oixy in the “horizontal” plane,
see text. These active ions are shown as small black spheres.
TABLE I: Correlation induced corrections to the diagonal
Hamiltonian matrix elements for the valence-band O 2s, 2p
and conduction-band Mg 3s, 3p states. All numbers are in
eV. Negative corrections induce upwards shifts of the valence
bands and shifts to lower energies for the conduction bands.
∆Hnn(0) O 2s O 2p Mg 3s Mg 3p
On-site orb. relaxation −2.64 −2.04 — —
NN orb. relaxation −1.23 −1.20 −0.81 −0.84
Long-range polarization −1.80 −1.80 −2.25 −2.25
Total −5.67 −5.04 −3.06 −3.09
wards shift of the center of gravity of the 3s−3p band
complex.
The data listed in Table I indicate that the on-site
orbital relaxation and relaxation and polarization effects
at the nearest oxygen sites in the presence of an extra
electron or extra electron hole results in a reduction of the
HF band gap by about 4.05 eV, that is, more than 45% of
the difference between the HF and experimental values.
Large corrections are also expected to arise from long-
range polarization effects. The long-range polarization
energy of a dielectric due to the presence of an extra
charge±e can be expressed as ∆E(∞) = ∆E(R)−C/R1,
where
ǫ0 − 1
e2 , (11)
ǫ0 is the static dielectric constant of the material, and R
defines a sphere around the extra charge beyond which
the dielectric response reaches its asymptotic value ǫ0.
The energy increment ∆E(R) denotes the relaxation and
polarization energy up to the radius R around the added
particle. The constantC can be obtained by choosing two
different radii R1 and R2 where the quantities ∆E(R1)
and ∆E(R2) are computed, see for example Ref. 35.
However, we adopt here a simpler approach. We cal-
culate the corrections due to long-range polarization by
using the experimental value for the static dielectric con-
stant, ǫ0 = 9.7. Since relaxation and polarization effects
related to the nearest oxygen neighbors were already ac-
counted for, both, for the valence-band hole states and
the conduction-band electrons (see Table I) and since the
core-like electrons of the Mg2+ ions can be ignored in
these calculations, we set R as the average of the radii of
the first and second oxygen coordination shells around a
localized 2p hole or 3s(3p) electron: R = (a
2/2 + a)/2
for the O 2p valence-band states, where a = 4.217 Å is
the lattice constant50, and R = (a/2 + a
3/2)/2 for the
conduction-band states. The corrections to the diagonal
matrix elements of the Hamiltonian are then
∆Hnn(0) = −
ǫ0 − 1
, (12)
about −1.80 eV for the O 2s/2p bands and −2.25 eV for
the Mg 3s/3p bands. These numbers are also included in
Table I.
Before we discuss band narrowing (or broadening) due
to correlation effects which also affect the band gap, we
consider the loss of ground-state correlations. The latter
leads again to a shift of the center of gravity of the bands.
As pointed out in the previous section, some of the config-
urations that are present in the N -particle ground-state
are blocked when an electron is added or removed. We
investigated such correlations by CI calculations with sin-
gle and double excitations (CISD) and discuss first dif-
ferential effects for the N and (N−1) states. Since the
oxygen valence-band Wannier orbitals are rather local-
ized, we designed a cluster with a single O ion in the
active region. Around this central O site we added one
shell of Mg ions (6 Mg’s) and two shells of anions (12+6
O’s) to build the buffer region CB. In the CISD calcula-
tions for the N and (N−1) configurations we correlate
the 2s and 2p orbitals of the central O ion. Thereby
the occupancy of the hole orbital is kept frozen in the
calculations for the (N−1) states, which is referred to
as the frozen hole approximation38,39. Sets of separately
optimized orbitals were used for the hole states of the
(N−1)-particle system, as discussed above. We found
that for a 2p hole the correction to the on-site matrix
element of the Hamiltonian is ∆Hnn(0) = 0.85 eV, i.e.,
the 2p valence bands are downshifted by 0.85 eV. For the
O 2s hole states, this correction amounts to 0.99 eV. One
would expect similarly an upwards shift of the conduction
bands. However, the situation is somewhat different here.
When an extra electron is attached to the Mg2+ ion, it
polarizes the closed shells of the core. This is the domi-
nant effect now because Mg2+ has no valence electrons.
We may employ for our analysis the high-quality results
obtained for a free Mg ion by Doll et al.53. The correc-
tion to the ROHF Mg+→Mg2+ ionization potential was
found to be 0.27 eV in Ref. 53. A similar differential cor-
relation effect is occuring for the conduction-band states
in bulk MgO. Therefore the conduction bands are shifted
downwards instead of upwards and a partial cancellation
between loss of ground-state correlation effects for the
valence and conduction bands is taking place. The net
result is a slight increase of the gap between the valence
and conduction bands, in the range of 0.5 eV. Among
the different contributions discussed here, this appears
to produce the smallest corrections to the gap. More
advanced calculations for studying such differential cor-
relation effects are left for future work.
To summarize the results listed in Table I, relaxation
and polarization effects in bulk MgO are responsible for
a reduction of the Hartree-Fock gap by 8.1 eV, which
represents about 95% of the difference between the HF
and experimental values, 16.2 and 7.8, respectively. Im-
proved agreement is expected between our results and
the experimental data when applying higher-quality ba-
sis sets. We mention in this context that a reduction of
3.8 eV is obtained for the HF gap when going from va-
lence double-zeta to triple-zeta basis sets. Nevertheless,
this large energy difference is mainly related to the very
poor representation of the conduction-band states in the
TABLE II: Nearest–neighbor (NN), RNN = (1, 1, 0)a/2, and
next-nearest-neighbor (NNN), RNNN = (1, 0, 0)a, hopping
matrix elements for the conduction-band Mg 3s and 3p or-
bitals, see text. Results of frozen-orbital CI (FO-CI) are listed
in the second column; NOCI results in terms of separately
optimized, relaxed orbitals (RO-NOCI) are given in the third
column. All numbers are in eV.
tnn′(R) FO-CI RO-NOCI
tNN :
3s − 3s 0.41 0.42
3px(y) − 3px(y) 0.66 0.69
3px(y) − 3py(x) 0.72 0.77
3pz − 3pz 0.13 0.13
tNNN :
3s − 3s 0.36 0.37
3px − 3px 0.77 0.74
3py(z) − 3py(z) 0.13 0.12
calculations with double-zeta basis functions and such ef-
fects will be less substantial by further extension of the
basis sets. We also keep in mind that differential correla-
tion effects due to the existence of a different number of
electrons in the system’s ground-state and in the (N±1)
excited states determine a small correction in the oppo-
site direction, i.e., a slight increase of the fundamental
B. Off-diagonal matrix elements
We discuss next the effect of correlations on the widths
of the different bands. For that purpose the off-diagonal
matrix elements of the effective Hamiltonian (10) have to
be determined, i.e., the so-called hopping terms. At the
Hartree-Fock level, these matrix elements are obtained
by solving 2×2 secular equations where both wavefunc-
tions are expressed in terms of localized HF orbitals, see
Eqs. (4), (5), and (8). Relaxation and polarization ef-
fects in the nearby surroundings of the added electron
(or hole) are obtained by separate SCF optimizations
for the (N±1) states. The separate optimization of the
(N±1) wavefunctions leads to sets of non-orthogonal or-
bitals. There will be thus both Hamiltonian and over-
lap matrix elements between the extra electron (extra
hole) wavefunctions ΨN±1
0nσ and Ψ
. The calculation
of such matrix elements has been recently implemented
in molpro54. It is based on the transformation of the
CI vectors to bi-orthogonal orbitals and follows an idea
suggested by Malmqvist55. A similar approach was devel-
oped by Broer et al. in Groningen56. Instead of putting
these matrix elements directly into the eigenvalue equa-
tions that determine the Bloch energies ǫn(k) we extract
from these data effective hopping parameters associated
with various pairs of orbitals. Comparison between such
effective hopping terms and the HF off-diagonal Hamil-
tonian matrix elements offers an insightful picture of how
correlation effects modify the inter-site interactions and
FIG. 5: Sketch of the [Mg28O36] cluster employed for the
calculation of the NN conduction-band hopppings. The active
region CA includes two NN Mg sites, small grey spheres, plus
the bridging and apical ligands, large black spheres.
consequently the widths of the bands. For energetically
degenerate states, the effective hopping is defined as
tnn′ = (Hnn′ − Snn′Hnn)/(1− S2nn′) , (13)
where Hnn′ and Snn′ are the Hamiltonian and over-
lap matrix elements between the (N ± 1) states n and
n′. Since the separately optimized wavefunctions are ex-
pressed in terms of sets of non-orthogonal orbitals, this
type of secular problem is usually referred to as non-
orthogonal CI (NOCI). In the case of mutually orthogo-
nal states Snn′ is zero and tnn′(R) = Hnn′(R).
Nearest-neighbor and next-nearest-neighbor (NNN)
hopping matrix elements for the more diffuse conduction-
band Mg 3s and 3p orbitals are listed in Table II. In the
rocksalt structure, there are two O ions bridging two NN
cations, whereas two NNN Mg’s share a single oxygen.
We designed a [Mg28O36] cluster for calculating the NN
hoppings and a [Mg36O47] cluster for the NNN matrix
elements. As active regions, we employed [Mg2O6] and
[Mg2O9] kernels, respectively, see Fig.5 and Fig.6. As
in the calculations for the on-site matrix elements, for
each of these clusters we included in the buffer region CB
all metal and O ions in the first two coordination shells
of the active oxygens. Since the second cluster consists
of 83 atoms, we reduced the computational effort by re-
moving the polarization functions at the oxygen sites for
this particular cluster. Two different values are given in
Table II for each hopping matrix element. Results ex-
tracted from 2×2 CI calculations in terms of frozen HF
Wannier orbitals are listed in the second column. In the
third column, we allowed for full relaxation of the 2s and
2p orbitals of the O ions included in the active region CA,
for each of the (N+1) configuration state functions en-
tering the 2×2 CI. As illustrated in Figs. 5 and 6, those
active anions are the ligand(s) bridging two Mg sites and
the ligands which are nearest oxygen neighbors of one of
the active Mg sites and also of the bridging ligand(s).
The results show that the separate optimization of the
(N+1) wavefunctions induces only minor changes on the
electron hoppings, in the range of few percent. An inter-
esting feature is that variations occur in both directions,
i.e., some of the effective hoppings are slightly enlarged
by taking into account short-range relaxation and polar-
ization effects and some are reduced. Since these changes
are quite small, the width (and the overall structure) of
the lower conduction bands will change very little. Re-
garding the longer-range polarization effects, we expect
that their influence on the hopping terms is negligible.
In the case of non-equivalent orbitals, neighboring Mg
3s and Mg 3p, it is more difficult to define some NOCI ef-
fective hoppings and the corresponding data is not shown
in Table II. Nevertheless, the effect of short-range relax-
ation and polarization is also small for these interactions,
with reductions of the CI splittings of few meV.
The same type of analysis was applied to the (N−1) O
2s1/2p5 hole states. The spatial extent of the oxygen or-
bitals and the inter-site matrix elements are significantly
smaller. Since the largest relaxation effects concern or-
bitals at the same site, see Table I, we included in a first
step only two ligands in the active regions of our clusters.
Nevertheless, few other atomic shells were added around
these active sites in the buffer region CB. We considered
a [O28Mg10] cluster for computing the relaxation effects
on the NN hoppings and a [O32Mg11] cluster for the NNN
terms. The results are collected in Table III. The largest
corrections arise for the NN matrix elements, with abso-
lute values that are similar to the corrections obtained for
the NN conduction-band hoppings. In relative numbers,
these corrections are somewhat larger, by 10% to 20%
for the NN hoppings. A very different situation occurs in
strongly correlated oxides such as the layered cuprates,
where the existence of an antiferromagnetic spin back-
ground determines a reduction of the effective quasipar-
ticle hoppings by a factor of four57.
FIG. 6: Sketch of the [Mg36O47] cluster employed for the cal-
culation of the NNN conduction-band hopppings. The active
region CA includes two Mg sites, small grey spheres, and nine
O neighbors, large black spheres.
TABLE III: NN, RNN = (1, 1, 0)a/2, and NNN, RNNN =
(1, 0, 0)a, hopping matrix elements for the valence-band
states. Results of frozen-orbital CI (FO-CI) are listed in
the second column; NOCI results in terms of separately op-
timized, relaxed orbitals (RO-NOCI) are given in the third
column. All numbers are in eV. Data for the oxygen 2s or-
bitals are also included, although the O 2s band is much be-
low the O 2p bands. The values in parentheses include orbital
relaxation effects at four additional O sites, see text.
tnn′(R) FO-CI RO-NOCI
tNN :
2s− 2s 0.10 0.12 (0.11)
2px(y) − 2px(y) 0.32 0.37 (0.36)
2px(y) − 2py(x) 0.42 0.49 (0.47)
2pz − 2pz 0.12 0.14 (0.13)
tNNN :
2s− 2s 0.01 0.01
2px − 2px 0.06 0.06
The trends displayed in Table III are confirmed by
multi-reference CISD [MRCI(SD)] calculations. We cor-
related in these calculations the eight 2s and 2p orbitals
at the two NN oxygen sites and the reference active space
included those two orbitals involved in the hopping pro-
cess. The hopping matrix element is half of the en-
ergy separation between the lowest two eigenstates. The
MRCI values for the 2px−2px and 2px−2py hoppings, for
example, are 0.33 and 0.43 eV. These numbers are again
larger than the HF results, although the magnitude of the
effect is smaller as compared to the NOCI calculations.
Having the data of the SCF calculations for the (N−1)
oxygen hole states at hand, we performed an analysis of
the changes produced in the composition of the (relaxed)
orbitals in the immediate vicinity of an O− 2s1 or 2p5
anion. We found that an oxygen hole is causing polariza-
tion and “bending” of the 2p orbitals at the NN ligand
sites. The “bending” of the NN 2p orbitals towards the
O− site takes place through both p−s and pi−pj mixing.
These effects result in stronger inter-site orbital overlap
and explain the fact that the effective NOCI hoppings
are larger than the corresponding HF values. A similar
analysis for the conduction-band (N+1) states is compli-
cated by the presence of several sets of “active” orbitals,
Mg 3s, 3p and bridging O 2p, and we could draw no clear
conclusions in that case.
Extra calculations were performed for the NN hop-
ping matrix elements with four additional oxygen sites
included in the active region of the cluster, CA. Those
are the four ligands which are nearest neighbors to both
O ions involved in the hopping process. They are situated
in the median plane of the segment RNN = (1, 1, 0)a/2.
A cluster composed of 62 sites, [O40Mg22], was employed
for these calculations. The results are given in Table III
in parentheses. The corrections due to relaxation and po-
larization at the four nearest O neighbors are very small,
0.01 to 0.02 eV. The effect of these corrections is to re-
duce somewhat the absolute values of the NN hoppings.
Experimental studies for characterizing the valence
electronic structure of MgO have been carried out
using x-ray photoelectron spectroscopy (XPS)58 and
angle-resolved ultraviolet photoelectron spectroscopy
(ARUPS)59. The measured width of the O 2p bands
is about 6.5 eV58,59. The HF valence-band width is 5.50
eV and inclusion of local correlations leads to a slight
broadening of the O 2p bands, which brings the ab initio
result in good agreement with the experiment. For com-
parison, density-functional calculations within the local
density approximation and using the same Gaussian ba-
sis sets as in the HF calculations predict a width of 4.68
eV for the O 2p bands.
The correlation induced corrections to the widths of
the bands also influence the band gap. In the fcc lat-
tice the dispersion of p bands at the Γ point depends on
two of the nearest-neighbor hoppings, ǫx(Γ) = const. +
(110)
x,x − 4t(011)x,x + ...60,61. With the notations from Ta-
ble III, t
(011)
x,x = t
(110)
z,z . Corrections of 0.04 eV for t
(110)
and 0.01 eV for t
(011)
x,x (or t
(110)
z,z ), see Table III, imply an
upwards shift of the O 2p bands at the Γ point and a
narrowing of the fundamental gap by about 0.30 eV. For
the Mg 3s−3p conduction-band complex such changes at
the Γ point are smaller because the correlation induced
corrections to the 3s−3s and 3s−3p inter-site matrix
elements are lower.
IV. SUMMARY AND CONCLUSIONS
We have analyzed the different correlation contribu-
tions to the energy gap of MgO and to the widths of
the conduction and valence bands. This was done within
the quasiparticle description. As regards correlation ef-
fects we have distinguished between relaxation and po-
larization around an electron or hole added to the (neu-
tral) ground-state. The net result is a reduction of the
Hartree-Fock gap from 16.2 eV to a value of 8.1 eV. This
has to be compared with a measured energy gap of 7.8
eV. Within the local density approximation (LDA) to
density functional theory a gap of 5.0 eV is found. This
is not surprising since LDA is known to produce too small
gaps for insulators.
The calculations were performed with triple-zeta basis
sets. Since on the Hartree-Fock level the calculated gap
differs for double- and triple-zeta basis sets by 3.8 eV, one
may consider the good agreement with the experimental
gap as somewhat fortuitous. Triple-zeta basis sets are
known, however, to produce reliable results in quantum
chemistry and therefore a further extension of the basis
set should keep the corrections small. Presently we are
not able to work with larger basis sets.
It was shown that a large contribution to the correla-
tion induced corrections to the fundamental gap comes
from on-site and nearest-neighbor relaxation, i.e., from
the immediate neighborhood of the added particle. But
also the long-range part of the polarization generated by
the extra particle contributes substantially to the reduc-
tion of the HF gap. This long-range part can be treated
in a continuum approximation thereby using the known
dielectric constant of MgO. The so-called loss of ground-
state correlations makes a small contribution in MgO.
The reason is that the conduction-band Wannier orbitals
have predominant Mg 3s or 3p character. The added
electron will essentially go thus to a Mg site where it
polarizes the closed 1s2, 2s2, and 2p6 shells. This effect
reduces the gap and counteracts the loss of ground-state
correlations which occurs when an electron is removed
(hole state) and therefore can no longer contribute to the
correlations of the remaining ones. Finally, also changes
in the widths of the bands influence the energy gap.
One surprising effect which we found is an enhance-
ment of the width of the valence bands when local corre-
lations are taken into account. This is the opposite one
expects since correlations lead usually to band narrowing
instead of broadening. It results from a slight bending
of nearest-neighbor 2p orbitals towards a O− site when
we account for correlations. The bending of the nearest-
neighbor ligand orbitals increases the wavefunction over-
lap and hence facilitates the hopping.
The present results prove the usefulness of band cal-
culations based on quantum chemical techniques. They
allow for well controlled approximations and a transpar-
ent interpretation of the different microscopic processes
which determine the size of the energy gap and the widths
of the bands.
We thank Dr. E. Pahl and Prof. J. Fink for fruitful
discussions.
1 P. Fulde, Electron Correlations in Molecules and Solids
(Springer-Verlag, Berlin, 1995).
2 N. H. March, Electron Correlation in Molecules and Con-
densed Phases (Plenum Press, New York, 1996).
3 P. Fulde, Adv. Phys. 51, 909 (2002).
4 R. McWeeny, Methods of Molecular Quantum Mechanics
(Academic Press, London, 1996).
5 T. Helgaker, P. Jørgensen, and J. Olsen, Molecular
Electronic-Structure Theory (Wiley, Chichester, 2000).
6 S. Suhai, Phys. Rev. B 27, 3506 (1983).
7 C.-M. Liegener, J. Chem. Phys. 88, 6999 (1988).
8 J. J. Ladik, Phys. Rep. 313, 171 (1999).
9 J.-Q. Sun and R. J. Bartlett, J. Chem. Phys. 104, 8553
(1996).
10 S. Hirata and S. Iwata, J. Chem. Phys. 109, 4147 (1998).
11 P. Reinhardt and J.-P. Malrieu, J. Chem. Phys. 109, 7632
(1998).
12 P. Y. Ayala, K. N. Kudin, and G. E. Scuseria, J. Chem.
Phys. 115, 9698 (2001).
13 A. Abdurahman, A. Shukla, and M. Dolg, Phys. Rev. B
65, 115106 (2002).
14 C. Pisani, M. Busso, G. Capecchi, S. Casassa, R. Dovesi,
L. Maschio, C. Zicovich-Wilson, and M. Schütz, J. Chem.
Phys. 122, 094113 (2005).
15 S. Hirata and R. J. Bartlett, J. Chem. Phys. 112, 7339
(2000).
16 R. Pino and G. E. Scuseria, J. Chem. Phys. 121, 2553
(2004).
17 M. Albrecht and J. Igarashi, J. Phys. Soc. Jpn. 70, 1035
(2001).
18 M. Albrecht and P. Fulde, Phys. Stat. Sol. B 234, 313
(2002).
19 C. Buth, U. Birkenheuer, M. Albrecht, and P. Fulde, Phys.
Rev. B 72, 195107 (2005).
20 W. Förner, Int. J. Quantum Chem. 43, 221 (1992).
21 W. Förner, R. Knab, J. Č́ıžek, and J. Ladik, J. Chem.
Phys. 106, 10248 (1997).
22 K. Fink and V. Staemmler, J. Chem. Phys. 103, 2603
(1995).
23 S. Hirata, R. Podeszwa, M. Tobita, and R. J. Bartlett, J.
Chem. Phys. 120, 2581 (2004).
24 C. M. Zicovich-Wilson, R. Dovesi, and V. R. Saunders, J.
Chem. Phys. 115, 9708 (2001).
25 U. Birkenheuer and D. Izotov, Phys. Rev. B 71, 125116
(2005).
26 A. Shukla, M. Dolg, P. Fulde, and H. Stoll, Phys. Rev. B
57, 1471 (1998).
27 A. Shukla, M. Dolg, P. Fulde, and H. Stoll, Phys. Rev. B
60, 5211 (1999).
28 A. Abdurahman, A. Shukla, and M. Dolg, J. Chem. Phys.
112, 4801 (2000).
29 P. Sony and A. Shukla, Phys. Rev. B 73, 165106 (2006).
30 J. Rubio, A. Povill, J. P. Malrieu, and P. Reinhardt, J.
Chem. Phys. 107, 10044 (1997).
31 H. Stoll, Phys. Rev. B 46, 6700 (1992).
32 For a review, see B. Paulus, Phys. Rep. 428, 1 (2006).
33 S. Horsch, P. Horsch, and P. Fulde, Phys. Rev. B 29, 1870
(1984).
34 J. Gräfenstein, H. Stoll, and P. Fulde, Chem. Phys. Lett.
215, 611 (1993).
35 J. Gräfenstein, H. Stoll, and P. Fulde, Phys. Rev. B 55,
13588 (1997).
36 M. Albrecht, P. Fulde, and H. Stoll, Chem. Phys. Lett.
319, 355 (2000).
37 V. Bezugly and U. Birkenheuer, Chem. Phys. Lett. 399,
57 (2004).
38 U. Birkenheuer, P. Fulde, and H. Stoll, Theor. Chem. Acc.
116, 398 (2006).
39 E. Pahl and U. Birkenheuer, J. Chem. Phys. 124, 214101
(2006).
40 E. Pahl, U. Birkenheuer, H. Stoll, and P. Fulde, manuscript
in preparation.
41 V. R. Saunders, R. Dovesi, C. Roetti, M. Causà, N. M.
Harrison et al., crystal 2000, University of Torino, Italy.
42 U. Birkenheuer, C. Willnauer, M. von Arnim, W. Al-
sheimer, and D. Izotov, Scientific Report, Max-Planck-
Institut für Physik komplexer Systeme, Dresden, 2002.
43 U. Birkenheuer, manuscript in preparation.
44 C. Roetti, R. Dovesi, M. von Arnim, W. Alsheimer, and U.
Birkenheuer, the crystal-molpro interface, Max-Planck-
Institut für Physik komplexer Systeme, Dresden, Germany.
45 P. Pulay, Chem. Phys. Lett. 100, 151 (1983); S. Saebø and
P. Pulay, Chem. Phys. Lett. 113, 13 (1985).
46 H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby, M.
Schütz et al., molpro 2006, Cardiff University, United
Kingdom.
47 S. S. P. Parkin, C. Kaiser, A. Panchula, P. M. Rice, B.
Hughes, M. Samant, and S.-H. Yang, Nature Materials 3,
862 (2004); S. Yuasa, T. Nagahama, A. Fukushima, Y.
Suzuki, and K. Ando, Nature Materials 3, 868 (2004).
48 S. Choopun, R. D. Vispute, W. Yang, R. P. Sharma,
T. Venkatesan, and H. Shen, Appl. Phys. Lett. 80, 1529
(2002).
49 R. C. Whited, C. J. Flaten, and W. C. Walker, Solid State
Commun. 13, 1903 (1973).
50 S. Sasaki, K. Fujino, and Y. Takeuchi, Proc. Japan. Acad.
B 55, 43 (1979).
51 M. I. McCarthy and N. M. Harrison, Phys. Rev. B 49, 8574
(1994); R. Dovesi, C. Roetti, C. Freyria-Fava, E. Aprá, V.
R. Saunders, and N. M. Harrison, Philos. Trans. R. Soc.
London Ser. A 341, 203 (1992).
52 J. Pipek and P. G. Mezey, J. Chem. Phys. 90, 4916 (1989).
53 K. Doll, M. Dolg, P. Fulde, and H. Stoll, Phys. Rev. B 52,
4842 (1995).
54 A. O. Mitrushchenkov and H.-J. Werner, Mol. Phys., in
press (2007).
55 P.-Å. Malmqvist, Int. J. Quantum Chem. 30, 479 (1986).
56 R. Broer and W. C. Nieuwpoort, Chem. Phys. 54, 291
(1981); Theor. Chim. Acta 73, 405 (1988).
57 L. Hozoi, S. Nishimoto, and C. de Graaf, Phys. Rev. B, in
press (2007) [also available on cond-mat/0611720].
58 S. P. Kowalczyk, F. R. McFeely, L. Ley, V. T. Gritsyna,
and D. A. Shirley, Solid State Commun. 23, 161 (1977).
59 L. H. Tjeng, A. R. Vos, and G. A. Sawatzky, Surf. Sci.
235, 269 (1990).
60 J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
61 N. W. Ashcroft and N. D. Mermin, Solid State Physics
(Thomson Brooks/Cole, 1976).
http://arxiv.org/abs/cond-mat/0611720
