the Hartree-Fock method and in Density-Functional theory. In the following
we will see an historical approach of this kind: the Hartree method.
6.1
The Hartree Approximation
The idea of the Hartree method is to try to approximate the wave functions,
solution of the Schrödinger equation for N electrons, as a product of single-
electron wave functions, called atomic orbitals.
As we have seen, the best
possible approximation consists in applying the variational principle, by mini-
mizing the expectation value of the energy E = ψ|H|ψ for state |ψ .
The Hamiltonian of an atom having a nucleus with charge Z and N electrons
is
¯
h2
Zq2
q2
H = −
∇2 −
e +
e
(6.1)
2m
i
r
r
i
e
i
i
ij
ij
where the sum is over pairs of electrons i and j, i.e. each pair appears only
once. Alternatively:
N −1
N
=
(6.2)
ij
i=1 j=i+1
It is convenient to introduce one-electron and two-electrons operators:
¯
h2
Zq2
f
e
i
≡ −
∇2 −
(6.3)
2m
i
e
ri
q2
g
e
ij
≡
(6.4)
rij
With such notation, the Hamiltonian is written as
H =
fi +
gij
(6.5)
i
ij
47
provided that g act symmetrically on the two electrons (definitely true for the
Coulomb interaction).
6.2
Hartree Equations
Let us assume that the total wave function can be expressed as a product of
single-electron orbitals (assumed to be orthonormal):
ψ(1, 2, . . . , N ) = φ1(1)φ2(2) . . . φN (N )
(6.6)
φi(1)φj(1) dv1 = δij.
(6.7)
Variables 1, 2, .., mean position and spin variables for electrons 1,2,...;
dvi
means integration on coordinates and sum over spin components. Index i labels
instead the quantum numbers used to classify a given single-electron orbitals.
All orbitals must be different: Pauli’s exclusion principle tells us that we cannot
build a wave function for many electrons using twice the same single-electron
orbital. In practice, orbitals for the case of atoms are classified using the main
quantum number n, orbital angular momentum
and its projection m.
The expectation value of the energy is
ψ|H|ψ
=
φ∗1(1) . . . φ∗N (N )
fi +
gij φ1(1) . . . φN (N ) dv1 . . . dvN
i
ij
=
φ∗i(i)fiφi(i) dvi +
φ∗i(i)φ∗j(j)gijφi(i)φj(j) dvidvj
i
ij
=
φ∗i(1)f1φi(1) dv1 +
φ∗i(1)φ∗j(2)g12φi(1)φj(2) dv1dv2
i
ij
(6.8)
In the first step, we made use of orthonormality (6.7); in the second we just renamed dummy integration variables with the ”standard” choice 1 and 2.
Let us now apply the variational principle to formulation (4.10), with the constraints that all integrals
Ik =
φ∗k(1)φk(1) dv1
(6.9)
are constant, i.e. the normalization of each orbital function is preserved:
δ
ψ|H|ψ −
k Ik
= 0
(6.10)
k
where k are Lagrange multipliers, to be determined. Let us vary only the
orbital function φk. We find
δIk =
δφ∗k(1)φk(1) dv1 + c.c.
(6.11)
48
(the variations of all other normalization integers will be zero) and, using the
hermiticity of H as in Sec.4.1.1,
δ ψ|H|ψ
=
δφ∗k(1)f1φk(1) dv1 + c.c.
+
δφ∗k(1)φ∗j(2)g12φk(1)φj(2) dv1dv2 + c.c.
(6.12)
j=k
This result is obtained by noticing that the only involved terms of Eq.(6.8) are those with i = k or j = k, and that each pair is counted only once. For instance,
for 4 electrons the pairs are 12, 13, 14, 23, 24, 34; if I choose k = 3 the only
contributions come from 13, 23, 34, i.e.
j=k (since g is a symmetric operator,
the order of indices in a pair is irrelevant)
Thus the variational principle takes the form
δφ∗k(1) f1φk(1) +
φ∗j(2)g12φk(1)φj(2) dv2 − kφk(1) dv1 + c.c. = 0
j=k
i.e. the term between square brackets (and its complex conjugate) must vanish
so that the above equation is satisfied for any variation:
f1φk(1) +
φ∗j(2)g12φk(1)φj(2) dv2 = kφk(1)
(6.13)
j=k
These are the Hartree equations (k = 1, . . . , N ). It is useful to write down
explicitly the operators:
¯
h2
Zq2
q2
−
∇2
e φ
e
k (1) +
φ∗
φj(2) dv2 φk(1) = kφk(1)
2m
1φk (1) −
j (2)
e
r1
r
j=k
12
(6.14)
We remark that each of them looks like a Schrödinger equation in which in
addition to the Coulomb potential there is a Hartree potential:
q2
V
e
H (r1) =
ρk(2)
dv2
(6.15)
r12
where we have used
ρk(2) =
φ∗j(2)φj(2)
(6.16)
j=k
ρj is the charge density due to all electrons differing from the one for which we
are writing the equation.
6.2.1
Eigenvalues and Hartree energy
Let us multiply Hartree equation, Eq(6.13), by φ∗(1), integrate and sum over k
orbitals: we obtain
k =
φ∗k(1)f1φk(1)dv1 +
φ∗k(1)φk(1)g12φ∗j(2)φj(2)dv1dv2.
k
k
k j=k
(6.17)
49
Let us compare this expression with the energy for the many-electron system,
Eq.(6.8). The Coulomb repulsion energy is counted twice, since each < jk > pair is present twice in the sum. The energies thus given by the sum of eigenvalues of the Hartree equation, minus the Coulomb repulsive energy:
E =
k −
φ∗k(1)φk(1)g12φ∗j(2)φj(2)dv1dv2.
(6.18)
k
<jk>
6.3
Self-consistent potential
Eq.(6.15) represents the electrostatic potential at point r1 generated by a charge distribution ρk. This fact clarifies the meaning of the Hartree approximation.
Assuming that ψ is factorizable into a product, we have formally assumed that
the electrons are independent. This is of course not true at all: the electrons
are strongly interacting particles. The approximation is however not so bad if
the Coulomb repulsion between electrons is accounted for under the form of an
average field VH , containing the combined repulsion from all other electrons on
the electron that we are considering. Such effect adds to the Coulomb attraction
exerted by the nucleus, and partially screens it. The electrons behave as if they
were independent, but under a potential −Zq2e/r + VH(r) instead of −Zq2e/r of
the nucleus alone.
VH (r) is however not a “true” potential, since its definition depends upon
the charge density distributions of the electrons, that depend in turn upon the
solutions of our equations. The potential is thus not known a priori, but it is a
function of the solution. This type of equations is known as integro-differential
equations.
The equations can be solved in an iterative way, after an initial guess of the
orbitals is assumed. The procedure is as follows:
1. calculate the charge density (sum of the square moduli of the orbitals)
2. calculate the Hartree potential generated by such charge density (using
classical electrostatics)
3. solve the equations to get new orbitals.
The solution of the equations can be found using the methods presented in
Ch.2. The electron density is build by filling the orbitals in order of increasing energy (following Pauli’s principle) until all electrons are “placed”.
In general, the final orbitals will differ from the starting ones. The procedure
is then iterated – by using as new starting functions the final functions, or with
more sophisticated methods – until the final and starting orbitals are the same
(within some suitably defined numerical threshold). The resulting potential VH
is then consistent with the orbitals that generate it, and it is for this reason
called self-consistent field.
6.3.1
Self-consistent potential in atoms
For closed-shell atoms, a big simplification can be achieved: VH is a central
field, i.e. it depends only on the distance r1 between the electron and the
50
nucleus. Even in open-shell atoms, this can be imposed as an approximation,
by spherically averaging ρk. The simplification is considerable, since we know
a priori that the orbitals will be factorized as in Eq.(2.9). The angular part is given by spherical harmonics, labelled with quantum numbers
and m, while
the radial part is characterized by quantum numbers n and . Of course the
accidental degeneracy for different
is no longer present. It turns out that even
in open-shell atoms, this is an excellent approximation.
Let us consider the case of two-electron atoms. The Hartree equations,
Eq.(6.14), for orbital k = 1 reduces to
¯
h2
Zq2
q2
−
∇2
e φ
e
1(1) +
φ∗
φ2(2) dv2 φ1(1) = 1φ1(1) (6.19)
2m
1φ1(1) −
2(2)
e
r1
r12
For the ground state of He, we can assume that φ1 and φ2 have the same spher-
ically symmetric coordinate part, φ(r), and opposite spins: φ1 = φ(r)v+(σ),
φ2 = φ(r)v−(σ). Eq.(6.19) further simplifies to:
¯
h2
Zq2
q2
−
∇2
e φ(r
e
1) +
|φ(r2)|2 d3r2 φ(r1) = φ(r1)
(6.20)
2m
1φ(r1) −
e
r1
r12
6.4
Code: helium hf radial
Code helium hf radial.f901 (or helium hf radial.c2) solves Hartree equations for the ground state of He atom. helium hf radial is based on code
hydrogen radial and uses the same integration algorithm based on Numerov’s
method. The new part is the implementation of the method of self-consistent
field for finding the orbitals.
The calculation consists in solving the radial part of the Hartree equation
(6.20). The effective potential Vscf is the sum of the Coulomb potential of the nucleus, plus the (spherically symmetric) Hartree potential
Zq2
ρ(r2)
V
e
scf (r) = −
+ VH (r),
VH (r1) = q2
d3r2.
(6.21)
r
e
r12
We start from an initial estimate for VH (r), calculated in routine init pot (
(0)
V
(r) = 0, simply). With the ground state R(r) obtained from such potential,
H
we calculate in routine rho of r the charge density ρ(r) = |R(r)|2/4π (note that
ρ is here only the contribution of the other orbital, so half the total charge, and
remember the presence of the angular part!). Routine v of rho re-calculates
the new Hartree potential V out(r) by integration, using the Gauss theorem:
H
r
Q(s)
V out
H
(r) = V0 + q2e
ds,
Q(s) =
ρ(r)4πr2dr
(6.22)
r
s2
max
r<s
where Q(s) is the charge contained in the sphere of radius s; rmax is the out-
ermost grid point, such that the potential has already assumed the asymptotic
value V0 = q2e/rmax, valid at large r.
1http://www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/F90/helium hf radial.f90
2http://www.fisica.uniud.it/%7Egiannozz/Corsi/MQ/Software/C/helium hf radial.c
51
The Hartree potential is then re-introduced in the calculation not directly
but as a linear combination of the old and the new potential.
This is the
simplest technique to ensure that the self-consistent procedure converges. It
is not needed in this case, but in most cases it is: there is no guarantee that
re-inserting the new potential in input will lead to convergence. We can write
V in,new(r) = βV out
H
H
(r) + (1 − β)V in
H (r),
(6.23)
where 0 < β ≤ 1. The procedure is iterated (repeated) until convergence is
achieved. The latter is verified on the norm, i.e. the integral of the square, of
V out(r) − V in(r). square.
H
H
In output the code prints the eigenvalue 1 of Eq.6.13, plus various terms of the energy, with rather obvious meaning except the term Variational correction.
This is
δE =
(V in
H (r) − V out
H
(r))ρ(r)d3r
(6.24)
and it is useful to correct3 the value of the energy obtained by summing the eigenvalues as in Eq.(6.18), so that it is the same as the one obtained using
Eq.(6.8), where eigenvalues are not used. The two values of the energy are printed side by side.
Also noticeable is the ”virial check”: for a Coulomb
potential, the virial theorem states that T
= − V /2, where the two terms
are respectively the average values of the kinetic and the potential energy. It
can be demonstrated that the Hartree-Fock solution obeys the virial theorem.
6.4.1
Laboratory
• Observe the behavior of self-consistency, verify that the energy (but not
the single terms of it!) decreases monotonically.
• Compare the energy obtained with this and with other methods: per-
turbation theory with hydrogen-like wave functions (E = −5.5 Ry, Sect.
D.1), variational theory with effective Z (E = −5.695 Ry, Sect. D.3),
best numerical Hartree(-Fock) result (E = −5.72336 Ry, as found in the
literature), experimental result (E = −5.8074 Ry).
• Make a plot of orbitals (file wfc.out) for different n and . Note that
the orbitals and the corresponding eigenvalues become more and more
hydrogen-like for higher n. Can you explain why?
• If you do not know how to answer the previous question: make a plot
of Vscf (file pot.out) and compare its behavior with that of the −Zq2e/r
term. What do you notice?
• Plot the 1s orbital together with those calculated by hydrogen radial
for Hydrogen (Z = 1), He+ (Z = 2), and for a Z = 1.6875. See Sect. D.3
if you cannot make sense of the results.
3Eigenvalues are calculated using the input potential; the other terms are calculated using
output potential
52
Chapter 7
The Hartree-Fock
approximation
The Hartree method is useful as an introduction to the solution of many-particle
system and to the concepts of self-consistency and of the self-consistent-field,
but its importance is confined to the history of physics. In fact the Hartree
method is not just approximate: it is wrong, by construction, since its wave
function is not antisymmetric! A better approach, that correctly takes into
account the antisymmetric character of the the wave functions is the Hartree-
Fock approach. The price to pay is the presence in the equations of a non local,
and thus more complex, exchange potential.
7.1
Hartree-Fock method
Let us re-consider the Hartree wave function. The simple product:
ψ(1, 2, . . . , N ) = φ1(1)φ2(2) . . . φN (N )
(7.1)
does not satisfy the principle of indistinguishability, because it is not an eigen-
state of permutation operators. It is however possible to build an antisymmetric
solution by introducing the following Slater determinant:
φ1(1)
. . .
φ1(N )
1
.
.
.
ψ(1, . . . , N ) = √
(7.2)
N !
.
.
.
φN (1) . . . φN (N )
The exchange of two particles is equivalent to the exchange of two columns,
which induces, due to the known properties of determinants, a change of sign.
Note that if two rows are equal, the determinant is zero: all φi’s must be
different. This demonstrates Pauli’s exclusion principle. two (or more) identical
fermions cannot occupy the same state.
Note that the single-electron orbitals φi are assumed to be orthonormal:
φ∗i(1)φj(1)dv1 = δij
(7.3)
53
where the “integral” on dv1 means as usual “integration on coordinates, sum
over spin components” . We follow the same path of Sec. (6.2) used to derive Hartree equations, Eq.(6.13). Since a determinant for N electrons has N ! terms, we need a way to write matrix elements between determinants on a finite paper
surface. The following property, valid for any (symmetric) operator F and
determinantal functions ψ and ψ , is very useful:
φ∗
φ
1
1(1)
.
φ∗1(N )
1(1)
.
φ1(N)
ψ|F |ψ
=
.
.
.
F
.
.
.
dv1 . . . dvN
N !
φ∗ (1)
.
φ∗ (N )
φ (1)
.
φ (N )
N
N
N
N
φ1(1) . φ1(N)
=
φ∗1(1) . . . φ∗N (N )F
.
.
.
dv1 . . . dvN
(7.4)
φ (1)
.
φ (N )
N
N
(by expanding the first determinant, one gets N ! terms that, once integrated,
are identical). From the above property it is immediate (and boring) to obtain
the matrix elements for one- and two-electron operators:
ψ|
fi|ψ =
φ∗i(1)f1φi(1) dv1
(7.5)
i
i
(as in the Hartree approximation), and
ψ|
gij|ψ =
φ∗i(1)φ∗j(2)g12 [φi(1)φj(2) − φj(1)φi(2)] dv1dv2
(7.6)
ij
ij
The integrals implicitly include summation over spin components. If we assume
that g12 depends only upon coordinates (as in Coulomb interaction) and not
upon spin, the second term:
φ∗i(1)φ∗j(2)g12φj(1)φi(2) dv1dv2
(7.7)
is zero is i and j states are different (the spin parts are not affected by g12 and
they are orthogonal if relative to different spins).
It is convenient to move to a scheme in which the spin variables are not
explicitly included in the orbital index. Eq. (7.6) can then be written as ψ|
gij|ψ =
φ∗i(1)φ∗j(2)g12 [φi(1)φj(2) − δ(σi, σj)φj(1)φi(2)] dv1dv2
ij
ij
(7.8)
where σi is the spin of electron i, and:
δ(σi, σj) = 0 if σi = σj
=
1 if σi = σj
In summary:
ψ|H|ψ
=
φ∗i(1)f1φi(1) dv1
(7.9)
i
+
φ∗i(1)φ∗j(2)g12 [φi(1)φj(2) − δ(σi, σj)φj(1)φi(2)] dv1dv2
ij
54
Now that we have the expectation value of the energy, we can apply the vari-
ational principle. In principle we must impose normalization constraints such
that not only all φi stay normalized (as we did in the derivation of Hartree’s
equation) but also all pairs φi, φj with same spin are orthogonal, i.e., a (tri-
angular) matrix ij of Lagrange multipliers would be needed. It can be shown
however (details e.g. on Slater’s book, Quantum theory of matter) that it is
always possible to find a transformation to a solution in which the matrix of
Lagrange multipliers is diagonal. We assume that we re dealing with such a
case.
Let us omit the details of the derivation and move to the final Hartree-Fock
equations:
f1φk(1) +
φ∗j(2)g12 [φk(1)φj(2) − δ(σk, σj)φj(1)φk(2)] dv2 = kφk(1)
j
(7.10)
or, in explicit form,
¯
h2
Zq2
q2
−
∇2
e φ
e
k (1)
+
φ∗
[φj(2)φk(1)
(7.11)
2m
1φk (1) −
j (2)
e
r1
r
j
12
−
δ(σk, σj)φk(2)φj(1)] dv2 = kφk(1)
The energy of the system, Eq. 7.9, can be expressed, analogously to the Hartree case, via the sum of eigenvectors of Eq.(7.11), minus a term compensating the double counting of Coulomb repulsion and of the exchange energy:
E =
k −
φ∗k(1)φ∗j(2)g12 [φk(1)φj(2) − δ(σj, σk)φj(1)φk(2)] dv1dv2.
k
<jk>
(7.12)
Eq.(7.11) has normally an infinite number of solutions, of which only the lowest-energy N will be occupied by electrons, the rest playing the role of excited
states. The sum over index j runs only on occupied states.
Let us carefully observe the differences with respect to Hartree equations,