R. Dovesi et al.

: The Periodic Hartree-Fock Method 63

phys. stat. sol. (b) 217, 63 (2000)

Subject classification: 71.15.±m; 71.20.Ps; 71.55.Ht; 75.25.+z; 75.50.Ee; S9.11

The Periodic Hartree-Fock Method

and Its Implementation in the Crystal Code
R. Dovesi (a), R. Orlando (a), C. Roetti (a), C. Pisani (a),
and V.R. Saunders (b)

(a) Dipartimento di Chimica IFM, UniversitÁa di Torino, via Giuria 5, I-10125 Torino
(b) CLRC Laboratory, Daresbury, Warrington, Cheshire, WA4 4AD, U.K.

(Received August 10, 1999)

The present chapter discusses the Hartree-Fock (HF) method for periodic systems with reference
to its implementation in the Crystal program. The HF theory is shortly recalled in its Closed Shell
(CS), Unrestricted (UHF) and Restricted open shell (RHF) variants; its extension to periodic sys-
tems is illustrated. The general features of Crystal, the periodic ab initio linear combination of
atomic orbitals (LCAO) program, able to solve the CS, RHF and UHF, as well as Kohn-Sham
equations, are presented. Three examples illustrate the capabilities of the Crystal code and the
quality of the HF results in comparison with those obtained with the Local Density Approximation
using the same code and basis set: NiO in its ferro-magnetic and anti-ferromagnetic structure,
trapped electron holes in doped alkaline earth oxides, and F-centres in LiF.

1. Introduction
For three decades the Hartree-Fock (HF) method has been the most popular scheme
for the investigation of the electronic structure of atoms, molecules and clusters, either
as such or as a well defined starting point for more sophisticated techniques. Since the
beginning of the seventies, theoretical chemists have succeeded in developing power-
ful, general purpose, universally adopted HF-based ab initio computer programs [1].
The features of the most successful codes are essentially the same: the linearized HF-
Roothaan equations are solved by using a few localized functions per atom (usually
indicated as ``atomic orbitalsº, AOs) as a basis set; in turn, the AOs are expressed as
linear combinations of Gaussian type functions (GTF) with appropriate exponents and
``contractionº coefficients; after obtaining the molecular orbitals (MO), eigenvectors of
the Fock matrix, through a self consistent field (SCF) procedure, the correlation cor-
rection to the ground state wave function and properties, and the description of ex-
cited states are usually performed by means of configuration interaction (CI) or per-
turbation techniques. A posteriori, theoretical quantum chemistry appears to have
followed quite a linear development from the original formulation of the HF equa-
tions by Fock, Hartree and Slater [2, 3], their linearized expression by Roothaan and
Hall [4, 5], Boys' anticipation of present day techniques [6 to 8], to the full implemen-
tation of efficient computational tools. This has required a huge effort and ingenuity
as concerns the definition of standard and well assessed basis sets [9, 10], the devel-
opment of powerful analytic integration techniques [11 to 13], the efficient treatment
of hundreds of thousand configurations in CI calculations [14], the new ideas for the
64 R. Dovesi et al.

solution of the correlation problem [15]. In parallel, a set of semi-empirical electronic

structure theories have been developed which, starting from the HF equations, on the
one hand drastically reduce the computational effort by disregarding or approximating
most of the bielectronic integrals, on the other incorporate correlation effects by way
of largely intuitive approximations. A survey of these methods can be found in
Ref. [16].
Density Functional Theory [17, 18] (DFT), that had played for long time a minor
role in molecular quantum chemistry, has become more and more popular in the last
years, and DFT related computational schemes are now available in standard ab initio
molecular codes [19].
The situation in computational solid state physics and chemistry is quite different
from the one sketched above for molecular quantum chemistry.
A great variety of formal schemes, computational methods and techniques are
adopted for the determination of the electronic structure of crystalline compounds; the
various proposals (many of them are presented in this book) differ in many respects,
including the basis set (numerical, plane waves, localized, mixed), solution techniques
of the basic equations, selected Hamiltonian and all its ingredients. The relative merits
and limits of the various proposals are difficult to assess at the moment, as many of
the computer codes are not available to the scientific community and are often in
rapid evolution. The main reason for this situation is that computational solid state
physics is a relatively young science with respect to molecular quantum chemistry: ab
initio methods appeared in the late seventies, and the first general, portable, publicly
available code, Crystal88 [20], was distributed only ten years ago; in solid state com-
putational physics, it is just at the beginning that ``filteringº process that in molecular
quantum chemistry led, for example, to abandon Slater and lobe type orbitals [21],
owing to the difficulties encountered in the analytic evaluation of many mono- and bi-
electronic integrals in the former case [22], and because of the poor numerical accu-
racy obtained for integrals involving high angular quantum number AOs with the lat-
ter basis type.
There are, however, two features that are common to most periodic computer pro-
grams:
a) the use of DFT based Hamiltonians, the only exception being the Crystal [23, 24]
code, that can solve both the HF and Kohn-Sham (KS) equations;
b) the prevalent use of numerical rather than analytical techniques.
As a matter of fact, the HF approach has never been very popular among solid state
physicists, for many reasons; probably the most important two are:
±± HF performs poorly for the electron gas, the simplest of all periodic systems; it
has very often been assumed that this remains true (or possibly gets worse) for most
properties of real solids.
±± The non-local exchange term makes the HF equations more difficult to solve than
the KS ones, where the exchange±correlation potential mxc …r; ‰rŠ† is simply a multiplica-
tive operator, no matter how complicated its determination.
Early attempts at formulating and implementing HF computational schemes for peri-
odic systems had scant success, in spite of the high quality and generality of some of
them [25 to 30]. Only in recent years has it become possible to formulate a fair judge-
ment about the usefulness of the HF approach in solid state physics, since the advent of
The Periodic Hartree-Fock Method 65

powerful computers and of general purpose computational schemes [20, 23, 24, 31],
which have allowed us to assess its performance for a variety of systems [32].
In this chapter we will recall the general features of the HF scheme, both in its
formulation for closed and open shell systems, and its generalization to periodic com-
pounds when a local basis is used. Explicitly referring to the implementation in the
Crystal program, we will shortly recall some of the most critical points of this imple-
mentation and show that all the relevant interactions are evaluated analytically. This
requires a much bigger implementation effort, but ensures a very high numerical accu-
racy, that is necessary in many cases, for example, in chemical reactions [33, 34] and in
the study of magnetic interactions [35 to 37], that are of the order of a fraction of a
kJ/mol.
The chapter ends up with a few applications; all of them refer to cases in which
paramagnetic species are involved; the UHF results will be compared with DFT data
obtained with the same basis set.

2. Hartree-Fock Theory
The non-relativistic Born-Oppenheimer Hamiltonian operator for a finite cluster of N
nuclei and m electrons has the form

b ˆ ÿ1 P r2 ÿ P P ZA P P
m m N m iÿ1
H 2 i ‡ jri ÿ rj jÿ1 ‡ N ; …1†
iˆ1 iˆ1 Aˆ1 jri ÿ Aj iˆ1 jˆ1

where Zj denotes the atomic number of the j-th nucleus, and

P P ZA ZB
N Aÿ1
N ˆ …2†
Aˆ1 Bˆ1 jA ÿ Bj

is the nuclear repulsion energy. In the above equations atomic units are used. The Ha-
miltonian of Eq. (1) is parametrically dependent only on the nuclear coordinates, A. A
trial HF wave function is a single determinantal function constructed by antisymmetriz-
ing a Hartree product of m single particle spin orbitals,
m 
b Q
Y ˆA …wi …ri † c…i†† ; …3†
i

(w and c denote the orbital and spin functions, respectively; Ab is the antisymmetrizing
operator) and is obviously an approximate solution to Eq. (1). The HF electronic wave
function is variationally the best among these single determinantal functions.

2.1 Restricted closed shell (CS) theory for molecular systems

In the CS [5, 6] theory, the same molecular orbital (MO) is used to define both an a- and
a b-spin orbital; the wave function is then constructed from md …md ˆ m=2† doubly occu-
pied orbitals
m 
b Q …w …r2iÿ1 a† w …r2i b† :
d
Y CS ˆ A i i …4†
i

The expectation value of the Hamiltonian with respect to the trial HF wave function,
an approximation to the total electronic energy, takes the form
66 R. Dovesi et al.

ECS ˆ hY CS j H jY CS i
md „
P
ˆ2 dr fw*i …r† ‰ÿr2 =2 ‡ v…r†Š wi …r†g
i
md „
P
‡ dr dr0 ‰w*i …r† wi …r†Š ‰w*j …r0 † wj …r0 †Š=jr ÿ r0 j
ij
md „
P
ÿ 12 dr dr0 ‰w*i …r† wj …r†Š ‰w*j …r0 † wi …r0 †Š=jr ÿ r0 j : …5†
ij

The above total energy expression is to be minimized with respect to variations in the
orbitals, subject to the orthonormality constraint among the orbitals,
„
w*i …r† wj …r† dr ˆ dij : …6†
The resulting HF equation takes the form (fb stands for the Fock Hamiltonian)
„
fb wi …r†  ‰ÿ…1=2† r2 ‡ v…r† ‡ dr0 r…r0 †=jr ÿ r0 jŠ wi …r†
md „
P
ÿ dr0 ‰w*j …r0 † wi …r0 †Š=jr ÿ r0 jŠ wj …r†
j

ˆ ei wi …r† ; …7†
where the electron charge density r…r† is defined in terms of the wj …r† functions as
follows:
P
md
r…r† ˆ 2 jwj …r†j2 …8†
j

and the factor two is due to the double occupation of the orbitals. The first three terms
(kinetic, nuclear attraction, inter-electronic Coulomb repulsion) in Eq. (7) coincide with
the corresponding ones which appear in the KS equations of DFT. The last term repre-
sents the ``exactº non-local exchange operator, that in the KS equation is substituted by
the effective exchange±correlation potential, mxc …r; ‰rŠ†. Like the KS equations, the HF
equations must be solved through a Self-Consistent (SC) procedure, since both the
Coulomb and the exchange operators depend on the set of functions, wi …r†.
If we express the molecular orbitals, wi , as linear combinations of nb local (and real)
functions jm (AOs, the basis set),
P
nb
wi ˆ cmi jm : …9†
m

Eq. (7) transforms into the following matrix equation:

FC ˆ SCE ; …10†
where F and S are the Fock and the overlap matrices, C and E the eigenvectors and
eigenvalues. The expression of the various contributions to F is given below
Fmn ˆ Tmn ‡ Zmn ‡ Bmn ‡ Xmn ; …11†
where
„
Tmn ˆ ÿ 12 jm …r† r2r jn …r† dr ; …12†
The Periodic Hartree-Fock Method 67

…
P
N ZA
Zmn ˆ jm …r† j …r† dr ; …13†
Aˆ1 jr ÿ Aj n
P
nb
Bmn ˆ Plr …mn j lr† ; …14†
lr

P
nb
Xmn ˆ ÿ 12 Plr …ml j nr† ; …15†
lr

and
„„
…mn j lr† ˆ jm …r† jn …r† jr ÿ r0 jÿ1 jl …r0 † jr …r0 † dr dr0 …16†
P
Pmn ˆ 2 c*mi cni ; …17†
i

P
nb
ECS ˆ 12 Pmn …Fmn ‡ Tmn ‡ Zmn † : …18†
mn

2.2 Unrestricted (UHF) and restricted Half-closed-shell (RHF)

Hartree-Fock theories

For the description of systems containing unpaired electrons (such as atoms, molecules
with an odd number of electrons, radicals) a single determinant is not, in the most
general case, an appropriate wave function. In order to get the correct spin eigenfunc-
tion of these systems, it is necessary to choose a linear combination of Slater determi-
nants (whereas, in closed shell systems, a single determinant always gives the appropri-
ate spin eigenfunction). The UHF and RHF approximations permit, however, to
describe open shell systems while maintaining the simplicity of the single determinant
approximation.
In the RHF approximation [38, 39] the wave function is a single determinant con-
structed from md doubly occupied orbitals (the same MO being used to define both an
a- and a b-spin orbital), and a further ms singly occupied MOs whose spins are all
aligned in parallel and which constitute the half-closed-shell. We therefore have
m ˆ 2md ‡ ms .
In the UHF theory [40] the spin orbitals are taken to be of the form

w…x† ˆ wa …r† a…s† or wb …r† b…s† : …19†

The double occupancy constraint allows the RHF approach to obtain solutions that
are eigenfunctions of the spin operator, Sb2 , whereas UHF solutions are formed by a
mixture of spin states (both are eigenfunctions of Sbz ). The greater variational freedom
allows the UHF method to produce wave functions that are energetically more stable
than the corresponding RHF ones; another advantage of the UHF method is that it
allows solutions with locally negative spin density (i.e. anti-ferromagnetic systems), a
feature that RHF solutions cannot possess.
The molecular orbital diagram for the UHF and RHF cases is shown in Fig. 1. We
refer to the original papers quoted above, or to molecular quantum mechanics text-
books for the definition of the expressions equivalent to Eqs. (10) to (18) in the UHF
and RHF cases.
68 R. Dovesi et al.

Fig. 1. Molecular orbital diagram for the Restricted Open Shell (RHF, left) and the Unrestricted
(UHF, right) Hartree-Fock methods

3. Hartree-Fock Theory for Periodic Systems; the Crystal Program

At first sight, the reformulation of the HF-Roothaan equations (Eqs. (10) to (18))
sketched above, to make them suitable for crystalline systems is a relatively easy task,
and has been proposed more than thirty years ago by various authors [26 to 28]. Start-
ing from a set of localized functions in the unit reference cell of the crystal, jm …r†
…m ˆ 1; :::; p†, it is possible to build Bloch functions (BF),
P
fm …k; r† ˆ exp …ik
G† jm …r ÿ G† ; …20†
G

where the sum is extended to all G vectors of the direct lattice, while k is a general
vector in the first Brillouin zone (BZ). The resulting set of equations is in one-to-one
correspondence with Eqs. (10) to (18).

3.1 The formal equations for closed shell systems

In the basis of the Bloch functions, the Fock matrix becomes
P ikG
Fmn …k† ˆ e Fmn …G† ; …21†
G

where Fmn …G† is the matrix element of the Fock operator between the m-th AO located
in the zero cell and the n-th AO located in the G cell. The row index can be limited to
the zero cell for translational symmetry. Matrices represented in the Bloch basis (or in
``k spaceº) take a block diagonal form, as Bloch functions are bases for the Irreducible
Representations (IR) of the Translation Group (TG); each block has the dimensions of
the AO basis in the unit cell, nb ,
F…k† C…k† ˆ S…k† C…k† E…k† : …22†
In principle the above equation should be solved for all the (infinite) k points of the
Brillouin zone; fortunately, only a finite and usually small subset of these blocks, corre-
sponding to a suitable sampling of k points, needs to be diagonalized, because interpo-
The Periodic Hartree-Fock Method 69

lation techniques can be used for eigenvalues and eigenvectors throughout the first
Brillouin zone.
Fmn …G† ˆ Tmn …G† ‡ Zmn …G† ‡ Bmn …G† ‡ Xmn …G† ; …23†
„
Tmn …G† ˆ jm …r† r2r jn …r ÿ G† dr ; …24†
…
PN P ZA
Zmn …G† ˆ jm …r† j …r ÿ G† dr ; …25†
Aˆ1 M jr ÿ A ÿ Mj n
nb P
P P
Bmn …G† ˆ Plr …G0 † …m0nG j lMr…M ‡ G0 †† ; …26†
lr G0 M

P
nb P P
Xmn …G† ˆ ÿ12 Plr …G0 † …m0lM j nGr…M ‡ G0 †† ; …27†
lr G0 M

„ 0 P
nb
Plr …G0 † ˆ 2 dk eikG c*lj …k† crj …k† q…EF ÿ Ej …k†† ; …28†
j

P
nb P
ECS ˆ 12 Pmn …G† …Fmn …G† ‡ Tmn …G† ‡ Zmn …G†† ; …29†
mn G

where EF is the Fermi energy, the integration in Eq. (28) extends to the first Brillouin
zone, and cj …k† and Ej …k† are eigenvectors and eigenvalues of the F…k† matrix.
The above equations, as such, are useless, because the G, G0 and M summations
extend to the infinite set of translation vectors; a strategy must then be specified for the
treatment of the infinite Coulomb and exchange series, as well as for the substitution of
the integral that appears in Eq. (28) with a weighted sum extended to a finite set of k
points. An accurate and efficient solution to these problems has been implemented in
the Crystal code [23, 24, 31, 32]. Some of the features of Crystal will be summarized
in the next section. Before moving to these technical aspects, we want however to
spend a few more words on the periodic UHF method, as all the applications discussed
in the last section concern paramagnetic species.

3.2 The UHF equations

In the UHF case, two sets of matrix equations must be solved self-consistently, for a
and b electrons,
F a …k† Ca …k† ˆ S…k† Ca …k† Ea …k† ; …30†
b b b b
F …k† C …k† ˆ S…k† C …k† E …k† ; …31†
where the F a …k†, F b …k† and S…k† matrices are obtained by Fourier transform from the
corresponding ``direct spaceº equivalent quantities as in Eq. (21).
We now define the total density (P) and spin density (Pspin ) matrices,

P…G† ˆ Pa …G† ‡ Pb …G† ; …32†

spin a b
P …G† ˆ P …G† ÿ P …G† ; …33†
where Pa …G† and Pb …G† are obtained as in Eq. (28) by using the a and b eigenvectors
obtained from Eq. (30) and (31), respectively.
70 R. Dovesi et al.

The F a …G† and F b …G† matrices are defined as follows:

a spin
Fmn …G† ˆ Fmn …G† ÿ Xmn …G† ; …34†
b spin
Fmn …G† ˆ Fmn …G† ‡ Xmn …G† : …35†
The F…G† matrix is defined as in Eq. (23), where the total density matrix P…G† defined
spin
in Eq. (32) is used in the Coulomb and exchange terms; Xmn …G† is defined as Xmn …G†
spin
in Eq. (27), where however the spin density matrix P …G† is used instead of the total
density matrix P…G†.

4. The Implementation of the Periodic HF Equations

in the Crystal Program
Crystal is an ab initio periodic computer program implemented by the Group of Theo-
retical Chemistry of the University of Torino (Italy) in collaboration with the Computa-
tional Group at the Daresbury Laboratory (UK). The first version was distributed in
1988 [20], and for about ten years it has been the only public general-purpose ab initio
periodic computer program; three improved versions of Crystal were subsequently re-
leased; [23, 41, 42] (the last one is in distribution since the beginning of 1999). More
than 200 copies of both the '92 and '95 versions have been distributed over the world.

4.1 The general features

We summarize here the main features of the Crystal98 program:
1. It is a ``direct spaceº program, in the sense that all the relevant quantities (mono-
and bi-electronic integrals, overlap and Fock matrices, see Eqs. (23) to (29)), are com-
puted in the configuration space. Just before the diagonalization step the Fock matrix is

Fig. 2. Scheme of the implementation of the UHF method in Crystal

The Periodic Hartree-Fock Method 71

Fourier transformed to reciprocal space (Bloch functions basis; see Eq. (21)), then the
eigenvalues and eigenvectors of the F…k† matrices are combined to generate the ``direct
spaceº density matrix for the next SCF cycle (Eq. (28)).
2. In its most recent version, Crystal can solve the HF as well as the KS equations;
as regards the latter, the most popular local and non local functionals are available, as
well as hybrid schemes, such as the so called B3-LYP, which combines the HF exchange
term with the Becke [43] and Lee-Yang-Parr [44] functionals according to the formula
proposed by Becke [45]. As regards the former, the CS, RHF and UHF options are
available. Schemes are also available, that permit to correct the HF total energy by
estimating the correlation energy a posteriori, integrating a correlation-only functional
of the HF charge density. This latter scheme has been shown to provide accurate bind-
ing energies for a large family of compounds [32].
The general structures of the program is illustrated in Fig. 2, for the case of the UHF
path. The scheme shows that the most relevant quantities (Fock and density matrices,
eigenvalues and eigenvectors) are duplicated, so that it is possible to differentiate a
from b orbitals. The two sections of the code, corresponding to a and b electrons, are
independent until the Fermi energy calculation and Fock matrix reconstruction. It is
then possible to force the system into a state with a particular Sz value, by imposing the
desired ma ÿ mb value when the crystalline monoelectronic energy levels are populated
at each cycle of the SC step.
Many steps of the calculation are common to the HF and DFT options (for example
the treatment of the Coulomb series, which are evaluated analytically, as discussed be-
low); the main difference concerns, obviously, the exchange (and correlation) contribu-
tion to the Hamiltonian matrix and total energy: in the HF case the exchange bielectron-
ic integrals are evaluated analytically, and the exchange series is truncated after a
certain number of terms, as discussed below. In the DFT calculations [46], the ex-
change±correlation potential is expanded in an auxiliary basis set of GTF, with even
tempered exponents. At each SCF iteration the auxiliary basis set is fitted to the actual
analytic form of the exchange±correlation potential, which changes with the evolving
charge density. For numerical integration, the atomic partition method proposed by
Becke [47] has been adopted, combined with the Gauss-Legendre (radial) and Lebedev
(angular) quadratures.
3. It can treat systems periodic in 0 (molecules), 1 (polymers), 2 (slabs) and 3 (crystals)
directions with similar accuracy; this permits the evaluation of energy differences such as
bulk-minus-molecule (lattice energy of a molecular crystal), bulk-minus-slab (surface en-
ergy), bulk-minus-chain (inter-chain interactions) with high accuracy, as well as energy
differences between crystals with different cell size, shape, and number of atoms.
4. All the crystalline space, layer and rod groups are available, simply by specifying
the label of the group according to the international crystallographic notations.
5. Several geometrical options are available, which permit an easy manipulation of
the cell (creating defects, distorting the cell, building supercells, cutting slabs from the
bulk, extracting molecules from the bulk, and so on).

4.2 The treatment of the infinite Coulomb and exchange series

As anticipated, the problem with the formal periodic equations (22) to (29), is that they
contain three infinite summations (G, G0 and M), that extend to all lattice vectors. We
72 R. Dovesi et al.

refer to previous papers for an exhaustive presentation of the strategy adopted in Crys-
tal for the treatment of the Coulomb [31, 48] and exchange [24] series. Here we sum-
marize the main ideas:
±± The general strategy is to evaluate the integrals analytically, whenever possible. The
core of the integral package in the Crystal code derives from the GAUSSIAN70 [1]
and ATMOL [49] integral packages; both have been heavily modified and generalized
for many reasons, and new parts have been added. As regards the latter point, the
following new kinds of integrals are required for the construction of the Fock operator:
a) multipolar integrals of a product of two GTFs;
b) bipolar expansions of two interacting charge distributions represented by a pro-
duct of two AOs each;
c) interaction of a charge distribution rmn with the infinite array of point charges and
higher multipoles of a charge distribution. These integrals are evaluated analytically
using Ewald type techniques through a combination of recursion relations involving
Hermite polynomials and spherical harmonics, so that all the integrals can be obtained
from the generalized error functions, as for molecular integrals [31, 48].
Specific routines are available for periodicity in one, two and three dimensions (see
Refs. [31, 50]).
±± The overlap between two GTFs decays exponentially with distance. This property
can be exploited for the truncation of the G and G0 sums in the Coulomb term (see
Eq. (26) and (29)) and for two of the exchange summations (Eq. (28)). In a similar way
the overlap, kinetic, nuclear attraction integrals can be reduced to the consideration of
a finite number of neighbours of the reference cell. The screening criteria are applied
both to the contracted AOs, and to the single Gaussians of a contraction. An efficient
sorting of the G vectors and evaluation of the overlap between shell couples allow a
rapid selection of the integrals to be computed.
±± The density matrix elements decay with distance between the involved functions. In
a localized basis set, the elements of the density matrix of an insulator decay exponen-
tially with the distance between the two centres, the larger the gap the faster the decay.
Density matrix in direct space decays to zero with distance for conductors, too, but
much more slowly [24, 51].
This behaviour is exploited for an efficient truncation of the exchange series in the
Fock matrix (Eq. (27)) and total energy (Eq. (29)) expressions. The two above equa-
tions show that the M summation is limited by the exponential decay of the product
jm jl and for a similar reason the jG ÿ M ÿ G0 j distance cannot be too large. These two
conditions imply that the G and G0 vectors cannot be very different, although their
moduli could be large. However, the exponential decay of the density matrix permits us
to disregard integrals involving G and G0 vectors with large moduli. The exchange ser-
ies are then truncated according to these criteria. Long range effects are not taken into
account, as they are negligible. Exchange bielectronic integrals are evaluated as such or
through a bipolar expansion (see below), when the two charge distributions do not
overlap.
±± The multipolar expansion and the Ewald method for the Coulomb series. Let us
define a partition of the cell charge distribution, for example in terms of Mulliken
atomic charges rA,
PP
r…r† ˆ rA …r ÿ A ÿ M† ; …36†
A M
The Periodic Hartree-Fock Method 73

where
P PP
rA …r ÿ A† ˆ Plr …G0 † jl …r† jr …r ÿ G0 † ÿ ZA d…r ÿ A† …37†
l2A r G0

is translationally invariant and the last term is the nuclear charge contribution. The
sum of the electron±electron and electron±nuclei contributions to the Fock matrix,
Rmn …G† ˆ Bmn …G† ‡ Zmn …G† (Eq. (25) and (26)), can then be written as
P P „„
Rmn …G† ˆ rmn …r1 ÿ G† jr1 ÿr2 jÿ1 rA …r2 ÿ A ÿ M† dr1 dr2 …38†
A M

with
rmn …r1 ÿ G† ˆ jm …r1 † jn …r1 ÿ G† : …39†
If the charge distributions rA and rmn do not overlap for any M, Eq. (38) can be written
as
P P P 
Rmn …G† ˆ g`m …A† Fmn`m …G; A ‡ M† ; …40†
A `m M

so that Rmn is expressed in terms of the multipole moments of rA (N`m is a normaliza-

tion factor),
„
g`m …A† ˆ rA …r2 ÿ A† N`m X`m …r2 ÿ A† dr2 …41†
and of the field integrals
„
Fmn`m …G; A ‡ M† ˆ jm …r1 † jn …r1 ÿ G† X`m …r1 ÿ A ÿ M† jr1 ÿ A ÿ Mjÿ2`ÿ1 dr1 :
…42†
where X`m is a real solid spherical harmonic.
The infinite M summation can be performed by using the Ewald technique, which
reduces it to a finite number of new integrals, that can be evaluated analytically [31].
The main advantage of Eq. (38), with respect to Eq. (25) and (26), is that a large
number of terms (typically, more than 1000 for a packed system like MgO with a dou-
ble-z basis set) are summed up (Eq. (37)) and treated at the same time with a corre-
sponding saving factor in computational effort. If rmn and rA do overlap for some M
vectors close to the origin (let us call this set MB ), the same Ewald technique is used.
A correction is however applied, consisting in subtracting by direct summation the con-
tribution from the MB set in Eq. (40), and adding the corresponding contribution from
Eq. (25) and (26). The resulting finite (possibly small) number of bielectronic Coulomb
integrals and nuclear attraction integrals is then evaluated with standard techniques; it
is, however, to note that a non negligible fraction of integrals corresponding to the MB
set can be approximated by a bipolar expansion of both charge distributions, because it
may happen that, although rA is penetrating rmn , some of the contributions to rA are
not, and can then be approximated (see Ref. [24], pp. 47 to 52).
±± Point symmetry. A crucial point for the speed and accuracy of the calculation is the
complete exploitation of symmetry. Point symmetry is exploited at various levels in Crys-
tal, both in the direct (AO) and reciprocal space (Bloch functions) representations.
As concerns the former aspect, two points are worth noting:
±± The contributions to the Fock matrix Fmn …G† are evaluated only for the set of
indices mnG belonging to the irreducible set, i.e. to the minimal set necessary for the
74 R. Dovesi et al.

generation of the complete matrix through the use of the point symmetry operations.
The saving factor in CPU time and space requirements can be as large as h, the num-
ber of the symmetry operations in the point group. We remind here that translational
symmetry has already been used for limiting the first index m to the zero cell.
±± Disk space requirements for the storage of the bielectronic integrals can be re-
duced, and the Fock matrix reconstruction at each step of the SCF stage can be acceler-
ated by storing symmetrized sums of bielectronic integrals [24, 52].
Point symmetry is also exploited in those parts of the program which work in the
Bloch function representation.
±± In principle the Fock matrix must be diagonalized at all the k points of the first
BZ. The eigenvalues are then used for the determination of the Fermi energy, and the
eigenvectors for building the density matrix of the next cycle. It is, however, easily
shown that the Fock matrix can be diagonalized only in the so-called irreducible part of
the BZ (IBZ); the full set of eigenvectors can then be obtained by rotation of the
eigenvectors of the IBZ.
±± The set of operators that do not move a k point of the IBZ form the so-called
little co-group [53] of the k point, and can be used to block-diagonalize the Fock
matrix F…k†, by building symmetry adapted crystalline orbitals [54, 55]. This point is
of minor importance when small unit cell systems are considered, as most of the k
point have no symmetry at all; however, when large unit cell, highly symmetric sys-
tems are considered, the exploitation of this particular aspect of symmetry becomes
essential, as the order of the little co-group of the few k points to be considered is
quite high [56].

5. Examples of Applications
Applications of the HF-Crystalline Orbitals SCF method as implemented in the Crys-
tal program in many different areas have been performed by the present authors since
1980, when the first studies devoted to diamond and graphite were published [57, 58]; a
complete list of our publications can be found at the web site http://www.ch.unito.it/ifm/
teorica/crystal.html. Here, we concentrate on three applications involving uncoupled
electrons, namely a transition metal oxide (NiO), a trapped hole in CaO bulk where a
Li is substituted for a Ca atom, an F-centre (an electron trapped in an anion vacancy)
in LiF. In all cases the UHF option has been used; for comparison, the calculations
have been repeated at the spin polarized LDA (LSDA) level.

5.1 The electronic and magnetic properties of transition metal compounds:

the case of NiO

NiO is a prototype compound in the family of transition metal oxides, and has been the
subject of extensive theoretical investigations [59 to 64]. From an experimental point of
view, these materials are classified as antiferromagnetic (AFM) ionic insulators. The
band structures, obtained at the UHF and LSDA level (we refer to previous papers for
the basis set description and other technical details [63, 64]), are reported in Fig. 3; the
figure includes the valence bands, and three conducting bands: two of them are the
empty (at the HF level) eg levels, the third is an anti-bonding oxygen p band. The two
band structures differ in many respects:
The Periodic Hartree-Fock Method 75

Fig. 3. HF and LDA NiO band structure; ferromagnetic solution

Fig. 4. NiO valence band structure obtained with Hamiltonians containing various mixing of the
LDA and HF exchange potential. The three highest bands are empty, and correspond to the eg
empty states and to an antibonding p state. The LDA solution is metallic
76 R. Dovesi et al.

(i) the HF solution presents a large gap, whereas the LDA solution is metallic; this
behaviour is a direct consequence of the local character of the exchange term, as de-
monstrated by Fig. 4, where the results of various ``hybridº solutions, obtained by mix-
ing in various percentages the HF and LDA exchange operator in the Hamiltonian, are
reported.
(ii) The highest part of the valence band structure is mainly contributed by p oxygen
states in the HF solution, whereas d states are predominant in the LDA solution.
As regards bulk properties, such as the lattice parameter, they are in line with pre-
vious experience: UHF overestimates (4.26  A), whereas LSDA underestimates (4.10  A)
the experimental result (4.17 A).
The electronic structure can be described in terms of Mulliken charges, as shown in
Table 1. It turns out that the UHF solution is very ionic, with spin moments and net
charges very close to the formal ones (‡2 for Ni, eight electrons in the d shell), and a
small and negative Ni±O bond population, indicating that the interaction between the
two ions is essentially electrostatic, with a short range repulsion contribution. The ferro-
magnetic (FM) and AFM solutions are very similar, apart from the small spin polariza-
tion on the oxygen atom in the latter case; this is in qualitative agreement with the very
small energy differences between the two solutions, that must imply very small differ-
ences in the wave function. The LSDA situation is quite different: the net charges are
about 0.5 electrons smaller, and when a solution corresponding to two uncoupled elec-
trons is imposed (Sz ˆ 1), the spin polarization involves also the oxygen atom for about
0.5 electrons, whereas in the UHF case the oxygen spin polarization was much smaller
(0.08 electrons). The LSDA FM and AFM solutions are quite different. The bond po-
pulation, though small, is positive.
The very different spin polarization of the oxygen obtained with the two Hamilto-
nians in the FM state is also evident from the spin density map (Fig. 5, left), whereas
it disappears, for symmetry reasons, in the AFM solution (Fig. 5, right). As the super-
exchange mechanism is mainly related to this oxygen polarization (see, for example,
the discussion in Ref. [37]) it is clear why the LSDA Hamiltonian produces much
higher energy differences between the FM and AFM states than UHF. The behaviour

Ta b l e 1
Net charges (q), spin moments (m) and bond populations (qAB ) in the ferromagnetic and
anti-ferromagnetic states of NiO evaluated according to a Mulliken partition. In the sec-
ond row the number of electrons in the Ni d shell is reported. Units: electrons

HF LDA

FM AF FM AF

Ni ‡1.87 ‡1.86 ‡1.55 ‡1.50

q dNi 8.098 8.102 8.387 8.443
O ÿ1.87 ÿ1.86 ÿ1.55 ÿ1.50
Ni 1.923 1.913 1.556 1.185
m dNi 1.919 1.912 1.550 1.185
O 0.077 0 0.432 0
qAB Ni±O ±0.032 ±0.032 ‡0.010 ‡0.018
The Periodic Hartree-Fock Method 77

of this quantity (DE) as a function of the lattice parameter in the UHF case is given
in Fig. 6: it shows the rapid fall-down to zero (that is: the FM and AFM states have
the same energy) for large Ni±O distance as a consequence of the exponential reduc-
tion of the short range inter-ionic repulsion; this behaviour has been documented
experimentally in the case of the KMF3 perovskites (M ˆ Mn,Ni), in quite good
agreement with the UHF results [37]. The LSDA DE value are much larger and show
a different dependence on the Ni±O distance: it is ten times bigger at 3.80 
A, about

Fig. 5. HF and LDA spin density maps in the FM and AFM states; iso-density lines are drawn at
intervals of 0.01 (Bohr)ÿ3 ; continuous, dashed and dash-dotted lines correspond to positive, nega-
tive and zero values, respectively
78 R. Dovesi et al.

Fig. 6. Energy difference between the ferro-

magnetic and anti-ferromagnetic solutions for
NiO as a function of lattice parameter. The two
curves refer to HF energies (diamonds) and en-
ergies corrected with a correlation-only density
functional (circles)

20 times bigger at the experimental geo-

metry, 25 times bigger at 4.60 A. In the
cases where a direct comparison with ex-
perimental results is possible, as for the
above mentioned KMF3 perovskites, the
UHF super-exchange constant extracted
from DE turns out to be about one third
of the corresponding experimental quan-
tity [37].

5.2 Trapped electron holes in alkali doped alkaline earth oxides

Paramagnetic centres are formed when alkaline earth oxides, doped with alkali metal
ions, are exposed to ionizing radiation at low temperature [65, 66]. EPR spectra give
evidence of the formation of alkali metal trapped electron holes in oxides such as MgO,
CaO and SrO [65 to 68], that are usually denoted as [M]0 . The interest in these materi-
als is connected with two distinct applications: 1) their possible use as insulators in high
radiation environments, as they can suppress radiation damage; 2) their catalytic activ-
ity in the formation of higher order hydrocarbons from methane, for possible exploita-
tion in industrial processes [69, 70].
Electron Nuclear Double Resonance (ENDOR) spectroscopy has been applied to
[Li]0 - and [Na]0 -doped isostructural oxides MgO, CaO and SrO for the accurate meas-
urement of their hyperfine structure, which has provided useful information about the
spin density distribution and geometry of the defects [71, 72].
The electronic structure and properties of trapped electron hole centres can also be
investigated theoretically with the UHF method. If a defect is localized and its interac-
tions with the rest of the crystal are short-ranged, it can formally be represented by a
periodic model, such as that adopted by the Crystal program, provided that a suitable
supercell can be found [73 to 75], i.e. a multiple of the unit cell, that is large enough to
make the mutual interaction among neighbouring defects negligible, but still in the
range of a computationally accessible problem. As far as the condition of non-interact-
ing defects is satisfied, the method is fairly accurate and size consistent: the numerical
accuracy of the implemented algorithms is such that the calculated energy per supercell
of a crystal is precisely proportional to the supercell volume [74]. An analysis of the
dependence of the formation energy (Ef ) of [Li]0 centres in MgO, for example, on the
supercell size (see Table 1 in Ref. [76]) shows that big supercells may be necessary to
compute Ef very accurately, but a supercell of 32 atoms (S32 ) allows the calculation of
The Periodic Hartree-Fock Method 79

the defect formation energy within a few percent error and to account for the most
important relaxation effects.
The UHF calculation provides very similar electronic structures for all these systems.
The unpaired electron, that is formed after the substitution of an alkali metal ion for a
cation in alkaline earth oxides, could in principle be found in two different electronic
configurations:
±± fractions of the unpaired electron almost symmetrically delocalized on the oxygen
ions nearest neighbours of M‡ ;
±± the M‡ ion coupled to one of the oxygen ions, forming a well localized electron
hole.
UHF calculations lead to the localized electron hole solution unequivocally, in agree-
ment with the experimental indications, and the localized configuration is preferred to
the delocalized one by about 300 kJ/mol.
As examples, we refer to bulk [Li]0 and [Na]0 centres in CaO, calculated with the S64
supercell at the HF equilibrium geometry for CaO (lattice parameter: 4.83  A; expt.:
4.79 A). We refer to Refs. [76, 77] for technical details of the calculation. The band
structure is particularly useful for the characterization of the defect. Figure 7 shows the
valence and conduction bands of [Li]0 in CaO, that are essentially due to the p-type
oxygen orbitals. The levels corresponding to the p electrons of Oÿ are split from the
valence band originated by the p atomic orbitals of the O2ÿ ions: the populated levels,
i.e. the doubly degenerate px ±py and the a-pz levels (z along the ionic pair axis) move
to lower energy values (a-pz is stabilized for as much as ÿ608 kJ/mol from the bottom
of the O2ÿ p valence band), and a b-pz level appears in the large band gap, due to the
electron hole. This picture is supported by the atomic net charges of M‡ and Oÿ (O1 )
reported in Table 2 and obtained from Mulliken population analysis, which are very
close to 1 jej. Also the spin moment of O1 is very close to one, indicating that the
unpaired electron is almost completely localized and only very small fractions of spin
density are found on the neighbouring ions. A more detailed description of electron

Fig. 7. Band structure of Li-doped CaO obtained with the S64 supercell for a- and b-spin states
80 R. Dovesi et al.

Ta b l e 2
Net atomic charges (q) and spin moments (m) evaluated according to a Mulliken parti-
tion of the periodic densities. Units: electrons. Labels of the atoms as in Fig. 8

M O1 O2

q m q m q m

[Li]0 in CaO 0.980 0.000 ÿ1.047 0.985 ÿ1.932 0.002

[Na]0 in CaO 0.998 ÿ0.002 ÿ1.053 0.986 ÿ1.928 0.003

charge and spin density is given in Fig. 8 for a [Li]0 centre in CaO. The unpaired elec-
tron is clearly localized on O1 with some residual spin density spread in the region
around the defect, oscillating from positive to negative values.
It results from energy minimization that the M‡ and Oÿ monovalent ions move in
opposite directions and increase their mutual distance, R (Table 3), as they are less
attracted by each other than by their neighbouring bivalent O2ÿ and Ca2‡ ions. This is
in reasonable agreement with the experimental indications, despite the rough approxi-
mations in the model used to derive R from primary experimental data [66, 71, 72].
The theoretical method allows us to distinguish the individual displacements of M and
O1 from their lattice positions: the M‡ ions are between 3 and 4 times more mobile
than Oÿ. This increase in the M±O1 distance with respect to the perfect lattice inter-
ionic distance (4.83 
A) is remarkable when M ˆ Li. The higher mobility of the small
Li‡ ion allows a stronger electrostatic interaction with the surrounding ions and corre-
sponds to a higher relaxation energy, which is nearly three times as large as for Na.
The displacements of the neighbouring oxygen ions from their lattice positions are not
as large as for M, but the corresponding relaxation energies are far from negligible,
especially for M ˆ Na (31 kJ/mol), whereas it is less important for Li (5 kJ/mol).

Fig. 8. Total electron charge (left) and spin (middle) density maps in the (100) plane through the
Li‡ ±Oÿ pair. On the right, the spin density function along the vertical path dashed in the spin map.
The separation between contiguous iso-density curves in the maps corresponds to 0.01 e/Bohr3 in
the total charge density and 0.005 e/Bohr3 in the spin map. Symbols as in Fig. 5
The Periodic Hartree-Fock Method 81

Ta b l e 3
Relaxation of M and O1, and M±O1 distances (R) in  A. The experimental values (Rexp )
reported were obtained with different interpreting models

DzM DzO R Rexp

CaO : [Li]0 0.44 ÿ0.11 2.91 2.86±±3.05 a †, 2.83  0:19 b †

CaO : [Na]0 0.25 ÿ0.08 2.74 2.23 to 2.46 c †
a
† Ref. [72]; b † Ref. [66]; c † Ref. [71].

A straight comparison between calculated and experimental data is possible for the
parameters that determine the hyperfine structure of the spectra [72]: the hyperfine cou-
pling constants a and b and the nuclear quadrupole coupling constant P. These constants
appear in the formulation of the spin Hamiltonian in terms of which the ENDOR spec-
tra are interpreted.
The isotropic (a) and anisotropic (b) hyperfine coupling constants (in MHz) are ex-
pressed for each nucleus N as [78]
2m0
aN ˆ gbe gN bn raÿb …0† ; …43†
3h
 
2m 1
bN ˆ 0 gbe gN be T11 ÿ …T22 ‡ T33 † ; …44†
3h 2
where the spin density raÿb at N, the elements of the hyperfine coupling tensor T and
the electron g factor are the only terms which depend on the electronic structure of the
system (in our calculation we use the free electron g to approximate the true g factor,
which is usually an acceptable approximation for our purpose). The other multiplicative
factors in Eq. (43) and (44) are all tabulated constants [78, 79] (h, the Planck constant,
be and bn , the electronic and nuclear magnetons, m0 , the permeability of the vacuum
and gN , the nuclear g factor). T is a tensor of rank two, that is obtained as the quadru-
pole moment of the spin density at N. Its generic element has the form:
…  2 
N P P spin r dij ÿ 3ri rj
T ij ˆ PmnG jm …r† jn …r ÿ G† dr ; …45†
mn G r5

where the origin of the reference system is at nucleus N and ri denotes the i-th compo-
nent of r. Tii in equation (44) are the elements of T in diagonal form, following the
convention that T11 is the maximum module component.
The other important parameters in the spectrum are the nuclear quadrupole coupling
efg
constants PN, that measure the coupling of the electric field gradient (qN ) with the
nuclear quadrupole moment (QN ) at N,
efg
3 e2 qN QN
PN ˆ ; …46†
h 4I…2I ÿ 1†
where I is the total nuclear spin quantum number, e the electron charge magnitude, QN
efg
is measured experimentally and qN is obtained from a tensor of the same form as T
where the spin density matrix in (45) is replaced by the total density matrix.
The calculated and experimental values of a, b and P for the coupling of the un-
paired electron to the M nucleus in M-doped CaO are compared in Table 4. The best
82 R. Dovesi et al.

Ta b l e 4
Spin Hamiltonian hyperfine isotropic (a) and anisotropic (b) coupling constants and nu-
clear quadrupole coupling constant (P) at nucleus M in doped CaO. Comparison be-
tween calculated (HF) and experimental [72] values. Units: MHz

a b P

M calc. exp. calc. exp. calc. exp.

Li ÿ0.819 ÿ2.472 1.468 1.317 0.001 0.009

Na ÿ2.538 ÿ9.145 1.558 1.877 0.449 0.446

agreement between the calculated and experimental results is obtained for b, with a
percentage error within 15% of the experimental value. The same error range for b has
been found also for MgO and SrO [77].
The calculated value of PN is surprisingly good for CaO when N ˆ Na, though it
must be remarked that the variability of the values of Q reported by different authors
is rather large and the calculation of P is obviously affected by the same uncertainty
(data in Table 4 were obtained for QLi ˆ ÿ0:03 and QNa = 0.14 barn, as reported in the
International Tables of Physics and Chemistry [79], Vol. 9, p. 31). In the Li case, both
experimental and calculated P values are negligible.
As regards a, an accurate estimation is a more difficult task, since it requires a very
precise determination of spin density at nuclei where it can be nearly zero. This is
particularly true in this case, where the unpaired electron is very well localized on Oÿ
and very low spin densities are found anywhere else in the crystal (see the spin mo-
ment values in Table 2. Nevertheless, it has already been emphasized that Fig. 8 reveals
a well defined pattern for the residual spin density, with an alternation of positive and
negative zones, that is not random; it is rather a consequence of spin polarization [80]:
the electrons in the defect region with the same spin (a) as the unpaired electron, are
stabilized by exchange and super-exchange interactions and, for this reason, they inter-
act with it in a slightly more favourable way. Therefore, not only the absolute value of
a is important, but also its sign, as it specifies if a- or b-spin electrons prevail at the
nuclei. In Fig. 8 (right) the profile of spin density in CaO along the Li±Oÿ axis is
shown in detail. What determines the value of aLi is the depth of the small negative
peak (about 5  10ÿ4 jej/Bohr3 ) at Li, i.e. almost three orders of magnitude less than
spin density at O1 (Oÿ ) and one order of magnitude less than that at O2 (incidentally,
the twin peak at O2 implies pz character, like for O1 ). The agreement of the calculated
aN with the experimental data obtained from ENDOR spectra (Table 4) is only semi-
quantitative, though the sign and the order of magnitude are correct, and the calculated
aN values scale approximately with the same trend as the experimental data, as can also
be seen for MgO and SrO in Ref. [77]. From the computational point of view, particu-
lar caution is required to the level of convergence of the SCF cycle which is needed for
a spin density of the order of 10ÿ4 jej/Bohr3 to be determined with the necessary preci-
sion. Convergence of the order of 1  10ÿ5 Hartree on total energy, which is adequate
for the calculation of most observables, may be insufficient to evaluate a correctly, as is
shown in Table 5 for Li-doped CaO. Attempts at improving the basis sets (either the
valence or the core functions) and altering the geometry of the defect, as well as better
approximating the Coulomb and exchange series or using a denser net in reciprocal
The Periodic Hartree-Fock Method 83

Ta b l e 5
Spin Hamiltonian hyperfine isotropic (a) and anisotropic (b) coupling constants and nu-
clear quadrupole coupling constant (P) at Li in doped CaO as a function of the level of
convergence of the SCF cycle, when the precision in the calculation of total energy is
10ÿx Hartree. Units: MHz

x aLi bLi PLi

5 ÿ0.952 1.468 0.001

6 ÿ0.920 1.468 0.001
7 ÿ0.854 1.468 0.001
8 ÿ0.819 1.468 0.001
9 ÿ0.809 1.468 0.001
10 ÿ0.807 1.468 0.001
11 ÿ0.805 1.468 0.001
12 ÿ0.805 1.468 0.001

space, do not lead to any significant improvement in the evaluation of aN and the poor
agreement between calculation and experiment is likely to be considered as a limit of
the method: the spin densities involved are so low that the lack of electron correlation
and the spin contamination intrinsic in the UHF approximation (the UHF wave func-
tion is an eigenfunction of the Sbz operator and not of the total spin momentum Sb2 ) may
be important.
In order to have some indications about the importance of electron correlation in the
determination of a, a comparison with results obtained with DFT based Hamiltonians
would have been useful. However, we did not succeed in localizing the electron hole
onto one oxygen ion with DFT methods. In particular, the LSDA unpaired electron is
nearly completely delocalized over the six oxygens around the defect, in contrast with
the experimental evidence. This result seems to confirm the tendency of LSDA, already
observed for NiO (Fig. 5), to delocalize unpaired electrons in tightly bound states onto
neighbouring atoms. In the case of NiO, however, stoichiometry limits charge transfer
from each Ni to a single O ion (the overall delocalization is about 25%), whereas, in the
present case, six O ions are available and delocalization can proceed further.

5.3 F-centres in lithium fluoride

The capability of the method to reproduce the experimental values of b and P with a
quite satisfactory accuracy and, at the same time, the difficulty in computing a as accu-
rately as b and P have stimulated the investigation of other paramagnetic species,
which are characterized by a stronger Fermi contact, such as F-centres in alkali halides,
and for which very accurate ENDOR measurements [81] are available in the literature.
Alkali halides are similar to alkali earth oxides as concerns crystal symmetry (same
space group) and ionicity, and the F-centres offer the advantage of a simpler electronic
structure than the electron holes in doped oxides.
As an example, we consider F-centres in lithium fluoride. The model is a 64 atom
supercell crystal at the experimental lattice parameter, from which an F atom has
been removed. The F-centre consists of an unpaired electron localized inside the va-
cancy, as can be seen in Fig. 9 and from the Mulliken population analysis (Table 6),
but the picture differs from that of the trapped holes in alkaline earth oxides in some
respects: the electron wave function in the vacancy has spherical symmetry and is much
84 R. Dovesi et al.

more diffuse; relaxation of the nuclei surrounding the defect is negligible in this case
(the Li ions around the vacancy are displaced by only 0.03  A far from it, with an en-
ergy gain of just 1.7 kJ/mol). The unpaired electron in the vacancy is obviously more
weakly bound than the electron of Oÿ in the trapped hole centre. A bound a-spin level
appears in the valence±conduction band gap in this case, whereas the unpaired electron
state of [M]0 centres was stabilized below the O2ÿ p valence band. As in the previous
examples, we compared HF with LDA results, using Dirac-Slater exchange potential
[82] and the Perdew-Zunger parametrization of the Ceperley-Alder free electron gas
correlation [83]. A discrete bound (a) and the corresponding unbound (b) levels ap-
pear in the gap also in the LSDA band structure, which is qualitatively similar, but
differ from UHF quantitatively: the band gap is 936 kJ/mol (UHF: 2130 kJ/mol) and the
energy difference between the defect a and b levels is 156 kJ/mol (UHF: 981 kJ/mol).
This great difference in stability of the F-centre corresponds indeed to different elec-
tron charge and spin density distributions (Fig. 9). Mulliken analysis (Table 6) assigns
almost exactly one electron to the UHF vacancy, while there are only 0.887 electrons in
the LSDA vacancy. At the atoms around the defect, the spin moments calculated with
both methods (Table 6) are nearly three orders of magnitude smaller than in the va-
cancy. Nevertheless, the spin density profiles in Fig. 9 show that the two sharp sym-
metric side peaks, which correspond to the Li nuclei, are about as high as the large
middle peak centred at the vacancy. Since the heights of the Li peaks measure the spin
density at those nuclei and the Fermi contact interaction at Li depends on it, the corre-

Fig. 9. HF (top) and LDA (bottom) total electron charge (left) and spin (middle) density maps in
the (100) plane through the vacancy. On the right, the spin density function along the vertical path
dashed in the spin map. Symbols and units as in Fig. 8
The Periodic Hartree-Fock Method 85

Ta b l e 6
Net atomic charges (q) and spin moments (m) evaluated with the HF and LDA approx-
imations according to a Mulliken partition of the periodic densities. Units: electrons. La-
bels of the atoms as in Fig. 9

F-centre Li F

q m q m q m

HF ÿ1.002 1.085 0.978 0.003 ÿ0.975 ÿ0.009

LDA ÿ0.887 0.951 0.986 0.007 ÿ0.961 0.001

Ta b l e 7
Spin Hamiltonian hyperfine isotropic (a) and anisotropic (b) coupling constants at first
(Li) and second (F) nearest neighbours of the F-centre in LiF, as calculated in the HF
and LDA approximations and measured experimentally with ENDOR spectroscopy
[81]. Units: MHz

Li F

a b a b

HF 38.571 3.244 76.375 11.257

LDA 41.811 3.150 105.56 13.720
expt. 39.0 3.2 105.9 14.9

sponding a constant value is large and positive. The same happens for the F ions
around the vacancy. Table 7 shows that, in this case, where we deal with spin densities
at the nuclei of the order of 10ÿ2 jej=Bohrÿ3 (two orders of magnitudes larger than in
the trapped hole case) the agreement between the calculated and experimental values
is satisfactory not only for b, but also for a, especially when the first nearest neighbours
(Li) are involved. For farther neighbours of the vacancy the UHF error in a and b
tends to increase (see Table 7 for the F ions), whereas the LSDA approximation seems
to perform better. As a matter of fact, the tendency of LSDA to delocalize unpaired
electrons, which prevents from describing well localized states (3d electrons in NiO and
2p electrons in [M]0 ) correctly, is particularly suitable in the present case, where the
unpaired electron occupies a very diffuse 1s state and leads to results in better agree-
ment with experiment than UHF.

6. Conclusions
In this chapter we have recalled the fundamentals of the HF method and discussed how
it can be used for the study of periodic structures. For the latter purpose we have made
reference to a specific implementation, that embodied in the Crystal public code.
While this is not the most general possibility (in particular, the choice of using localized
Gaussian functions to describe the crystalline orbitals pervades the very structure of the
code and dictates many technical solutions), we have been able in this way to show
concretely some of the subtleties which are needed to realize from the basic equations
a powerful, efficient and useful tool.
Our second main objective has been to demonstrate through examples the usefulness
of the HF method in solid state applications. In this sense we have been necessarily
86 R. Dovesi et al.

confronted with DFT, the alternative approach which represents nowadays the favour-
ite choice for the overwhelming majority of workers in this field of studies. In spite of
all its merits, DFT is known to meet with serious difficulties in the description of sys-
tems with well localized spin densities which cannot be described through statistical
local averages of exchange interactions. The examples here reported (which have been
partly elaborated specially for the present work) are all concerned with this kind of
problems. In our opinion, they prove convincingly the merits of the HF technique if not
else as an auxiliary technique to DFT calculations.
There is, however, a more fundamental issue which justifies the effort towards the
continuous improvement of HF-based techniques, and this comes from molecular quan-
tum mechanics. At variance with DFT, standard quantum-chemical schemes can pro-
vide (in principle, at least) results of any required accuracy concerning ground and
excited states of many-electron systems. In this respect, programs like Crystal provide
a unique opportunity to implement and test post-HF schemes for crystalline systems.

Acknowledgements The present work is part of a project coordinated by A. Zecchina

and co-financed by the Italian MURST (Cofin98, Area 03). Financial support from the
Italian C.N.R. is acknowledged.

References
[1] D.J. Hehre, W.A. Lathan, M.D. Newton, R. Ditchfield, and J.A. Pople, GAUSSIAN70, Pro-
gram number 236, QCPE, Indiana University, Bloomington (Indiana) 1970.
[2] V. Fock, Z. Phys. 61, 126 (1930).
[3] J.C. Slater, Phys. Rev. 35, 210 (1930).
[4] C.C.J. Roothaan, Rev. Mod. Phys. 23, 69 (1951).
[5] G.G. Hall, Proc. Roy. Soc. A205, 541 (1951).
[6] S.F. Boys, Proc. Roy. Soc. A200, 542 (1950).
[7] S.F. Boys, Proc. Roy. Soc. A201, 125 (1950).
[8] S.F. Boys, G.B. Cook, C.M. Reeves, and I. Shavitt, Nature (London) 178, 1207 (1956).
[9] E. Clementi and R. Roetti, Atomic Data Nucl. Data Tables 14, 177 (1974).
[10] S. Huzinaga, Gaussian Basis Sets for Molecular Calculations, Elsevier Publ. Co., Amsterdam
1984.
[11] V. R. Saunders, in: Computational Techniques in Quantum Chemistry and Molecular Physics,
Eds. G.H.F. Diercksen, B.T. Sutcliffe, and A. Veillard, Reidel, Dordrecht 1975 (p. 347).
[12] L.E. McMurchie and E.R. Davidson, J. Comput. Phys. 26, 218 (1978).
[13] V. R. Saunders, in: Methods in Computational Molecular Physics, Eds. G.H.F. Diercksen and
S. Wilson, Reidel, Dordrecht 1983 (p. 1).
[14] V.R. Saunders and J.H. VanLenthe, Mol. Phys. 48, 923 (1983).
[15] M.F. Guest and S. Wilson, Chem. Phys. Lett. 75, 66 (1980).
[16] P. Deak, in: Properties of Crystalline Silicon, Ed. R. Hull, EMIS Datareview Series, Vol. 20.,
INSPEC, London 1999 (p. 245).
[17] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[18] W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 (1965).
[19] M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheese-
man, T.A. Keith, G.A. Peterson, J.A. Montgomery, K. Raghavachari, M.A. Al-Laham, V.G.
Zakrzewski, J.V. Ortiz, J.B. Foresman, J. Ciolowski, B.B. Stefanov, A. Nanayakkara, M. Chal-
lacombe, C.Y. Peng, P.Y. Ayala, W. Chen, M.W. Wong, J.L. Andres, E.S. Reploge, R. Gomperts,
R.L. Martin, D.J. Fox, J.S. Binkley, D.J. Defrees, J. Baker, J.P. Stewart, M. Head-Gordon,
C. Gonzalez, and J.A. Pople, GAUSSIAN94, revision C.3, Gaussian, Inc., Pittsburg (P.A.) 1995.
[20] R. Dovesi, C. Pisani, C. Roetti, M. CausaÁ, and V.R. Saunders, Crystal 88 Program number
577, QCPE, Indiana University, Bloomington (Indiana) 1989.
[21] H. Preuss, Z. Naturf. 11a, 823 (1956).
The Periodic Hartree-Fock Method 87

[22] C.A. Weatherford and H.W. Jones (Eds.), ETO Multicenter Molecular Integrals, Reidel, Dor-
drecht 1982.
[23] V. R. Saunders, R. Dovesi, C. Roetti, M. CausaÁ, N.M. Harrison, R. Orlando, and C.M. Zico-
vich-Wilson, Crystal98, User's Manual, UniversitaÁ di Torino, Torino 1998.
[24] C. Pisani, R. Dovesi, and C. Roetti, Hartree-Fock ab initio Treatment of Crystalline Systems,
Lecture Notes Chem., Vol. 48, Springer-Verlag, Berlin 1988.
[25] R. McWeeny, Proc. Phys. Soc. (London) 74, 385 (1959).
[26] J.M. AndreÂ, L. Gouverneur, and G. Leroy, Internat. J. Quantum Chem. 1, 451 (1967).
J.L. BreÂdas, J.M. AndreÂ, J.G. Fripiat, and J. Delhalle, Gazz. Chim. Ital. 108, 307 (1978).
[27] G. Del Re, J. Ladik, and G. BiczoÂ, Phys. Rev. 155, 997 (1967).
[28] A. Karpfen, Internat. J. Quantum Chem. 19, 1207 (1981).
[29] A.B. Kunz, Phys. Rev. B 6, 606 (1972).
[30] R.N. Euwema, D.L. Wilhite, and G.T. Surratt, Phys. Rev. B 7, 818 (1973).
[31] V.R. Saunders, C. Freyria-Fava, R. Dovesi, L. Salasco, and C. Roetti, Mol. Phys. 77, 629 (1992).
[32] R. Dovesi, in: Quantum-Mechanical ab initio Calculation of the Properties of Crystalline Ma-
terials, Ed. C. Pisani, Lecture Notes Chem., Vol. 67, Springer-Verlag, Berlin 1996 (p. 179).
[33] M. Catti, G. Valerio, R. Dovesi, and M. CausaÁ, Phys. Rev. B 49, 14179 (1994).
[34] M. Catti, F. Freyria Fava, C. Zicovich, and R. Dovesi, Phys. Chem. Minerals 26, 389 (1999).
[35] M.D. Towler, R. Dovesi, and V.R. Saunders, Phys. Rev. B 52, 10150 (1995).
[36] M. Catti and G. Sandrone, Faraday Discuss. 106, 189 (1997).
[37] R. Dovesi, F. Freyria Fava, C. Roetti, and V.R. Saunders, Faraday Discuss. 106, 173 (1997).
[38] R. McWeeny, Proc. Roy. Soc. (London) A241, 239 (1957).
[39] I.H. Hillier and V.R. Saunders, Internat. J. Quantum Chem. 4, 503 (1970).
[40] J.A Pople and R.K. Nesbet, J. Chem. Phys. 22, 571 (1954).
[41] R. Dovesi, V. R. Saunders, and C. Roetti, Crystal92, User's Manual, UniversitaÁ di Torino,
Torino 1993.
[42] R. Dovesi, V. R. Saunders, C. Roetti, M. CausaÁ, N.M. Harrison, R. Orlando, and E. ApraÁ,
Crystal95, User's Manual, UniversitaÁ di Torino, Torino 1996.
[43] A.D. Becke, Phys. Rev. A 38, 3098 (1988).
[44] C. Lee, W. Yang, and R.G. Parr, Phys. Rev. B 37, 785 (1988).
[45] A.D. Becke, J. Chem. Phys. 98, 5648 (1993).
[46] M.D. Towler, A. Zupan, and M. CausaÁ, Comp. Phys. Commun. 98, 181 (1996).
[47] A.D. Becke, J. Chem. Phys. 88, 2547 (1988).
[48] R. Dovesi, C. Pisani, C. Roetti, and V.R. Saunders, Phys. Rev. B 28, 5781 (1983).
[49] D. Moncrieff and V.R. Saunders, Doc NAT648, University of Manchester Regional Comput-
ing Center, 1986.
[50] V.R. Saunders, C. Freyria-Fava, R. Dovesi, and C. Roetti, Comp. Phys. Commun. 84, 156
(1993).
[51] M. CausaÁ, R. Dovesi, R. Orlando, C. Pisani, and V.R. Saunders, J. Phys. Chem. 92, 909
(1988).
[52] R. Dovesi, Internat. J. Quantum Chem. 29, 1755 (1986).
[53] M. Lax, Symmetry Principles in Solid State and Molecular Physics, John Wiley & Sons, New
York 1974.
[54] C.M. Zicovich-Wilson and R. Dovesi, Internat. J. Quantum Chem. 67, 299 (1998).
[55] C.M. Zicovich-Wilson and R. Dovesi, Internat. J. Quantum Chem. 67, 311 (1998).
[56] C.M. Zicovich-Wilson and R. Dovesi, Chem. Phys. Lett. 277, 227 (1997).
[57] R. Dovesi, C. Pisani, F. Ricca, and C. Roetti, Phys. Rev. B 22, 5963 (1980).
[58] R. Dovesi, C. Pisani, and C. Roetti, Internat. J. Quantum Chem. 17, 517 (1980).
[59] T. Oguchi, K. Terakura, and A.R. Williams, Phys. Rev. B 28, 6443 (1983).
[60] K. Terakura, T. Oguchi, A.R. Williams, and J. KuÈbler, Phys. Rev. B 30, 4734 (1984).
[61] A. Svane and O. Gunnarson, Phys. Rev. Lett. 65, 1148 (1990).
[62] S. Massida, M. Posternak, and A. Baldereschi, Phys. Rev. B 46, 11705 (1992).
[63] W.C. Mackrodt, N.M. Harrison, V.R. Saunders, N.L.Allan, M.D. Towler, E. ApraÁ, and R. Do-
vesi, Phil. Mag. 68, 653 (1993).
[64] M.D. Towler, N.L. Allan, N.M. Harrison, V.R. Saunders, W.C. Mackrodt, and E. ApraÁ, Phys.
Rev. B 50, 5041 (1994).
[65] G. Rius, R. Cox, R. Picard, and C. Santier, C.R. Acad. Sci. (France) 271, 824 (1970).
88 R. Dovesi et al.: The Periodic Hartree-Fock Method

[66] O.F. Schirmer, J. Phys. Chem. Solids 32, 499 (1971).

[67] H.T. Tohver, B. Henderson, Y. Chen, and M.M. Abraham, Phys. Rev. B 5, 3276 (1972).
[68] M.M. Abraham, Y. Chen, J.L. Kolopus, and H.T. Tohver, Phys. Rev. B 5, 4945 (1972).
[69] D.J. Driscoll and J.H. Lunsford, J. Phys. Chem. 89, 4415 (1985).
[70] T. Ito, J. Wang, C.H. Lin, and J.H. Lunsford, J. Amer. Chem. Soc. 107, 5062 (1985).
[71] G. Rius and A. HerveÂ, Solid State Commun. 15, 399 (1974).
[72] M.M. Abraham, W.P. Unruh, and Y. Chen, Phys. Rev. B 10, 3540 (1974).
[73] C. Freyria-Fava, R. Dovesi, V.R. Saunders, M. Leslie, and C. Roetti, J. Phys.: Condensed
Matter 5, 4793 (1993).
[74] R. Orlando, R. Dovesi, P. Azavant, N.M. Harrison, and V.R. Saunders, J. Phys.: Condensed
Matter 6, 8573 (1994).
[75] R. Orlando, P. Azavant, M.D. Towler, R. Dovesi, and C. Roetti, J. Phys.: Condensed Matter
8, 1123 (1996).
[76] A. Lichanot, C. Larrieu, R. Orlando, and R. Dovesi, J. Phys. Chem. Solids 59, 7 (1998).
[77] A. Lichanot, C. Larrieu, C. Zicovich-Wilson, C. Roetti, R. Orlando, and R. Dovesi, J. Phys.
Chem. Solids 59, 1119 (1998); Erratum: J. Phys. Chem. Solids 60, 855 (1999).
[78] J.A. Weil, J.R. Bolton, and E. Wertz, Electron Paramagnetic Resonance ±± Elementary Theo-
ry and Practical Applications, John Wiley & Sons, New York 1994.
[79] D.R. Lide (Ed.), Handbook of Chemistry and Physics, 72nd ed., CRC Press, Boca Raton
1991/1992.
[80] O.F. Schirmer, J. Phys. C 6, 300 (1973).
[81] B. Henderson, Defects in Crystalline Solids, Edward Arnold, Ltd., London 1972.
[82] P.A.M. Dirac, Proc. Cambridge Phil. Soc. 26, 376 (1930).
[83] J.P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).